


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1994 


Evaluation of energy-sink stability criteria for 
dual-spin spacecraft 


Ortiz, Vincent Michael. 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/28185 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NPS community. 
get Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS'‘s first 
KNOX appointed — and published — scholarly author. 





LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library Monterey, California USA 93943 


















ae) ee | . 

















tie tesa Ape - a! 
. Ou EdD ee ‘ eos 
er P 7 ‘ s¢ “|. 
rr are ees Hag Bee 
ogre e 9 v 'e ’ " i) 
: oa 30s 6 ee lve the? 
eRe a ae ‘ 
i anes Soe ae he as 
sore 2 . itare iis nh as 
Lise Hat tac ' . to Ar Lt * M ; i. 
a ; ; re ae Phan nee apes Roce ‘ rte wie 
' wou no uA Pi Le Klar te [ee tome mas It, "Aer vey 
‘ 1 H ae of tat avails ¢. 
; Up ve ' ‘ . . aad U 








Ye CP Sl ay Ts 
Ag Gast clan 

‘ He oe ietslios 

Seep me ik Te a ae 


‘ f 
a an % Th way 0 wht pis 1,3 t 















+ dee 
' 


























































































































































Sane PT ere bidet itite Ole gt tee! tae 
: ‘oe umadee, te of Part mun Pry fs 
De Shirin sia ° Vd beget ce cages 4d thet fa 
ianghetantits wal ii Pt ts if, Soaps 
. ‘ ama, 4 4, eo a Rhee dhe 
Pee Y on tate. ‘gif 4 .'seteal? £4 #F aden 
“oe Mt AES Cid, gitede ter ant PANS lag. 
| ' 4 < + dary ase 1 
dad tae o ge it ma tt Z°Ghe'f Motshig h toate ith. fy ut 
‘ tei Me : Con ty LST Oe ests Po ty ‘ 8 eft Sos 
; 3 ST tee thpege: eFtes 
A . % ota uu ah Steet a Halts te 
ji ‘ ; 4, ' 3 hinge a 
a mes Heian to 
. i] Haeptes sy! pes 'y Vito ope Re : Halwa reeta s Ry apes | 
4 * * . s ‘ ‘ at 
Aer vit ifs 4 tie nent oea"b font see an “ ate 
‘ taera t ‘ es edates ay ry 
As idl es tan fy 
rr te 
ot 
bs | 
i] 
‘ slar she 
ae toe ae Ast Rea iatob ged 
. ees ‘ad be ellenet hel yi 
: Sou ‘a ' + 
orig POR Lar Cio f A ee | a 
ee | Pare Weal 
‘ re eT ; a pyr 
rer Van” ee TN | is . ov 
1, : ; 
au . ' 1 Pee PY x Died dgtatet ne 
F oe . ‘ eat at Ne Pb Seales cas re 
Vue ie tee ' repeat ar el Vo ee “4 vue sybetae 
' $ '* 3 2 
t . ‘ 5 ! ’ a _ ee . a4 
a : ot ¢ 4 Here, 4 : arias ate : ae 1. ‘ede atcn nate 
H 60; Pie ' : a Den LMC beer i tare fret) at m4 aoe Fatal g heh: 
. Ci {eae eT Reso et : Fee eh ea LO era rT any = H ead arth: Piers Oe 
1 tee ay ardcg ‘ a rer tac s u nie ' or ia Ra Ye Fe ne Fame W165 Rete. 
on ' ioe see LAUREL aia De is nes ‘ we Veddscactes 
. . 6 oe i ate ae et d =a i Om | a. Cre ia J f ¢ acact 
' Fs 1 ' " Aen le ‘ De his 
: " . ' PCr wy ieee aie ake a neared Gone e 
. oon a A .? a ' oP eo. A et Le ae * wt, 









os 
=% 5 
wha 





ts° 


7 D ‘ 
or vee Mate the 


and 











~@F ar 































F r 
ee ie aah avd 4, 
: 1 2 
F bee 
F H 
ae TUS ata ‘ 
= 0 dn4ud tin ° Pt ee 
te bet My), nhac ON Se tee 
1" id 
Tey gtee Een 
i mrt whi, 
. t, eye," 
t + 





Bore ' 
@ o.t gae 
oe fet) 5, 










thetav . 
t: Wedise 1 Se 
a de ween gy 



















. nates 
4 orate S.Somg’ 
tr. 


= 





ee” ew 
. fo 







fees 
‘ Lhe mug 

opto 
4 vg?! 















a” 





-o" 

Sages 
‘ Pen 
we. on Se te 
sao pa 
Sate ty e Bs 














vrs 
de 


3 
“ 


“a8 
od FU Lek 





* 
. 
by - 
bd . ? 
. ay . . st 2 8 te 
— neeee ey 4 ' ' Ceo eT BAC) . ae = ee ‘ 
ae ie eee aie) bn s ty De caglt > Aes Palen Rc eat peach Ue 
' ' ' pe PS “4 . e 4 ae bs Ta : ") tae %, 
‘ ’ ee a? Mo Bay a sohe bd wo 68 4g 
ae =f JF Bye to c ‘ 1 tes ' rT : ‘ © hy 
n ’ ° Lie Cale "of «6 ’ ep 
Sate : earn ad . ai? ty oy ’ Cie hata ser 
= . . a . . ese ' , . 
: ’ os wt * ary ‘ Less ae ’ eas cy aie 
id ‘ . o. e . aNe Pha) + 4 ae 
- ou . ‘ ‘ ee 
on 8 y 2 ie C . oe -4¢ 
-,! et ‘ + 4 * 
Coe . 













8 Por) ‘a4 












































































































































































































































































2% Cheat 2 + 
isla’ © op eae Garr, 
: * s oh re 
‘ 
$548 Loe 4 7; « rr 
a“ some g Ut 4g a7 Meee i gegidse ak finesse 
1 CF i 4 
' PUTA Leal ae ae 
; oe NOR es taka Rae 
' . . in Sexsee ° 
. . ‘ e i 
J ' . ' at Fee teen t 
a S "se 8 teeee zi 
‘ 2. ft glengaaes 
1 1 4 
. . 
Ls See 
: ies . ' 
. ' . ) £ gt 
' + a 
a aE oly 
. ethane 
’ m 3 ‘ “5 yt ati, 28 , 
' anit on 734 140 
re ha TS el 
; a tae of Siw irs ot ar | st yendosss, 
ae Rebs ra aL orfereiat ae etre 
: terstag fed Sp. ee Sto’ af eabets oo 
. Fame y Sn 1d S88 HS “df ee woe pte ss 
! eee SEEM, ee cital pie done ogee ot 
eg aS <. et SATU ty Met ate® Pa acliucge At atcha Y 
7 ! ‘ : Wn pet ro et rs ee Ti zayh, $ate ete A per gene 4: 
' s . tee Pee tee at * 
. A a's 4 fay ; 4 
' “e nhele een ace wc) 
: ' “tse 
. oe ea. 1 "es, 8 
' i at . i anf aug 
, cary aa acre ertataite « 
» ?¢ ' ' he ate ee ty sae 
ote asl, Bie tetera a Re 
i ' oe {3728 ara 
. ' a3 ere due Te a 3 wt 
aoe *. 48 .? 
_ =< U fea 
: aut . 
ae goatee a ed 
' 
: *. 
1 . s ar *4 ae 
id 1f@an Pie nig 
’ : ‘ "4 Cache pe rae 
‘ ; LT I) 
’ 
[A en | . 
’ "23 at tA 
. ’ « 
' 
: ' Chita etigy oke y 
Pa . ANS Coke @erbie gs Ns 6a ign 
7 ee keep Mier ea Ly ghee 
’ ray ele: TF yp "Dah ye ye, epg 
4 : nate platy ats Pare Lae Hi tot voter hgbarts edge se oe 
' 4 ' “fa sheets Oe ye PEAT S Bete ees ey 
de Og atan , = Pade iad ety bf a 
- Ps: ind PALS area ete 
_— eet rata ly? 
. Po ‘a 
. oe . \y 
‘ 
. | Bay 
“ te ey | 
' “€ ont peta 6 EM 
F ; ? Bera ha ieee wag ‘es Fah inital 
Ady g 4 8 Poy 1% 
4989 rg wlasy peak 
a : : stats Nall ec khey 
+f 2 t oe Ly ay ng 
: he ed ; ' a | % 
r ‘ Prete eetn, ees Pure rete e 
' ~ wg gt vey Hig) ost [eyed yt 
. oe n ree ‘ ah 4 “ae % K 
; ' ° , : 
‘ 4 a 2 i 
t 48 306 a os, 33 
as RUE an ae vente: 
‘ oe 5 es ; baad b, SG on? 5 
eed Af ia {isitaset48 
. ; UV Mot Fe Rat np oa Oy Bt iecrte 
aot Ay ee tas fp 
@- 4 ot * 
' f, § ty 
' 50; = 
' 
' 
' 










Pale 
wet eae 
$ sb pentahy shat: 
er "stp se, Pat 










eam tptzss frogs 
eee pe red te Ge eho] 


DUDLE a 
MONTE 33943-5101 








Approved for public release; distribution is unlimited 


Evaluation of Energy-Sink Stability Criteria 
For Dual-Spin Spacecraft 


by 
Vincent Michael va 
Lieutenant Commander, Uhtted States Navy 
B.ME., Georgia Institute of Technology, 1983 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN PHYSICS 
from the 


NAVAL POSTGRADUATE SCHOOL 
June 1994 


Unclassified 


Security Classification of this page 





| REPORT DOCUMENTATION PAGE 
la Report Security Classification UNCLASSIFIED 1b Restricuve Markings 





2a Securi Cinasification Authori 3 Distribution Availability of Report i 
2b _Declassification/Downgrading Schedule Approved for public release; distribution is unlimited 
4 Performing Organization Report Number(s) S$ Monitoring Organization Report Number(s) 

6a Name of Performing Organization 6b Office Symbol Ja Name of Monitoring Organization 

Naval Postgraduate School (If Applicable) RO Naval Postgraduate School 

6c Address (city, state, and ZIP code) 7b Address (city. state, and ZJP code) 

Monterey, CA 93943-5000 Monterey, CA 93943-5000 


8a Name of Funding/Sponsormg Orgenizatios | 8b Office Symbol 9 Procurement Instrument Identification Number 
(lf Applicable) 


8c Address (city, state, and ZIP code) 


10 Source of Funding Numbers 


11 Tobe (Include Security Classification) 
EVALUATION OF ENERGY-SINK STABILITY CRITERIA FOR DUAL-SPIN SPACECRAFT (U) 


12 Personal Author(s wa me ee 


incen 
13a Type of Report 13b Time Covered 14 Date of Report (year, moruh, day) 15 Page Count 
ht py Gs ; From To ne 5 94 152 
16 Supplementary Notation [all unclassified theses} | The views expressed m this thesis are those of the author and do not rellect 
the official policy or position of the Department of Defense or the U. S. Government 


eS a Se 
| | ae 


19 Abstract (continu on reverse if necessary and idenify by block asmber 



























18 Subject Terms (continue on reverse if necessary and identify by block aumber) 


Energy-Sink, Dual-Spin, Gyrostat 
Energy-Sink Stability Criteria 


The nutational stability of a dual-spin, quasi-rigid, axisymmetric spacecraft containing a driven rotor is analyzed. The 
‘purpose is to examine a revised energy-sink stability theory that properly accounts for the energy contribution of the motor. 
An inconsistency in the development disproves the existing energy-sink theory's assumption that the motor of the system 
contributes exactly enough energy to offset the frictional losses between the rotor and the platform. Using the concept of 
core energy, the revised stability criteria for a dual-spin, quasi-rigid, axisymmetnc spacecraft containing a driven rotor is 
denved. An expression for nutation angle as a function of core energy over time is then determined. Numerical simulations 
are used to verify the revised energy-sink stability theory. The dual-spin, quasi-rigid, axisymmetric system presented by D. L. 
Mingori was chosen for the simulation. Equations for angular momentum and total energy were necessary to validate the 
numerical simulation and confirm aspects of the revised energy-sink stability theory. These equations are derived from the 
first principles of dynamics and are included in the analysis. An explicit relationship for core energy as a function of time 
does not exist. Various models postulating core energy are presented and analyzed. The numerical simulations of the 
computed nutation angles as a function of the postulated core energy compare well with the actual nutation angles of the 
system to confirm the revised energy-sink stability criteria. 


21 Absoeci Secunty Classificauon 





Unclassified 
22a Name of Responsible Individual 22b Tel Include Area code 22c Office Symbol} 
I. Michael Ross (408) 36-2312 AA- 
DD PORM 1473, & MAR 83 APR edition may be used until exhausted 


All other editions are obsolete Unclassified 


ABSTRACT 
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Numerical simulations are used to verify the revised energy-sink stability theory. The dual- 
spin, quasi-rigid, axisymmetric system presented by D. L. Mingori was chosen for the 
simulation. Equations for angular momentum and total energy were necessary to validate 
the numerical simulation and confirm aspects of the revised energy-sink stability theory. 
These equations are derived from the first principles of dynamics and are included in the 
analysis. An explicit relationship for core energy as a function of time does not exist. 
Various models postulating core energy are presented and analyzed. The numerical 
simulations of the computed nutation angles as a function of the postulated core energy 
compare well with the actual nutation angles of the system to confirm the revised energy- 
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center of mass 
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moment vector of all forces acting on the system 

moment about the b; th coordinate axis 

mass of the rotor of the dual-spin system 
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total mass of the dual-spin system 
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respect to the reference frame B along the b2 axis 

rate of work done by the motor torque 
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distance the system center of mass is located along the b3 
axis from the rotor coordinate axes (b;) reference frame 
origin B* 

mathematical variable as defined in Equation (22) 

small perturbation value 

time rate of change of perturbation 

quantity defined in Equation (138) 

time rate of change of ¢€ 

postulated nutation angle as a function of rotor core energy 
postulated nutation angle as a function of platform core 
energy 


time rate of change of postulated nutation angle as a function 
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N oP 


Wo 


Oo 


Wj 


Q; 
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nutanion angle 
time rate of change of nutation angle 


rotor nutation frequency 

platform nutation frequency 

inertial nutation frequency 

quantity defined in Equation (138) 

quantity defined in Equation (138) 

quantity defined in Equation (138) 

relative rotation rate of the rotor with respect to the platform 
time rate of change of the relative rotation rate of the rotor 
with respect to the platform 

angular position of the platform with respect to the rotor, is 
equal to - of 


angular velocity vector of body and reference frame B with 
respect to inertial reference frame NV 

initial angular velocity value along the b3 th coordinate axis 
before the perturbation 

time rate of change of initial angular velocity value along the 
b3 th coordinate axis 

angular velocity of body B along the b; th coordinate axis 
time rate of change of the angular velocity of body B along 
the b; th coordinate axis 

second derivative with respect to time of the angular velocity 
of body B along the b; th coordinate axis 

spin angular velocity of body B, spin angular velocity of 
rotor of the dual-spin system 

time rate of change of spin angular velocity of body B, ttme 
rate of change of spin angular velocity of rotor of the dual- 
spin system 

spin angular velocity of platform of the dual-spin system 
time rate of change of spin angular velocity of platform of 
the dual-spin system 

transverse angular velocity vector 

transverse angular velocity vector magnitude 


xi 


W, time rate of change of transverse angular velocity vector 
magnitude 
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Mr = M+M’+4mp+ 4m,’ 


y a (a4’ + 4m,’) 
sa mani 
= 
t 4 _ m’ 
Poe vR 
_ ,(M+4mp) 
Lh Sage, M; 
(m’ + 4m’) 
Lo = L———_ 
2 Mz 
¢ = (pz +9'z’ 
A = 1, +1; + 2mpa2+2my/a” +(M’ + 4my/)(1 — VIL? 
C = 13+1]3'+4m,a2 + 4my’a” 


J3° = [3 + Bera OR 


A, = -(A-C)ano; 


Az = J3 0a 

A3 = 2MrlCa, 
2 

A4g = Mr¢ 


As = -Mrl’ ana; 
Ag = ~2m(C+Ly)z 


A7 = 2m(C+La)zana; 
Ag = ~2m(C+Lp)z0, 


xiv 


Ag = —2mzCa, 


Aio = 2mzzQ, 
Ai, = mz2 

Ai2 = —mz*@)Q; 
A13 = —MZQ 


A14 = —M2aQ Q? 

Ais = -2m'(C-Li}z’ 
Aig = 2m'(C- Ly)z' a0; 
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I. INTRODUCTION 
A chronological review of early spacecraft types, and the stability criteria developed 
for them, provides a sufficient background for the fundamental concepts of this thesis. The 
revised energy-sink stability criterion is then presented, and an equation for nutation angle 


as a function of energy dissipation over time is developed. 


A. SINGLE SPIN SATELLITES 
1. Equations of Motion 
The earliest satellites took advantage of the fact that stability could be achieved via 
the ‘gyroscopic stiffness’ of a spinning body. A preliminary dynamic model for the 
satellite could be achieved with Euler's equations of motion for a rigid body. Eu' 


moment equation can be written as 


_ “dh - 8dh MGB 
M cy ent ® xh (1) 


In component form this becomes 


M, = hy + @2h3- 3 h2 


h2 + @3 hh, — @, hy (2) 


M2 
M3 = h3 + @1h2-@2 hy 
Simplification of the model is accomplished by assuming that it is axisymmetric, that the 
body-fixed axes coincide with the principal axes (correct for simple spacecraft with 
I; = 12), and that the body is in a torque-free environment (a valid approximation used 


throughout this thesis). The spacecraft is then represented by Figure (1) and the equations 


of motion are 


0 = 1 + (3 -h) a o3 
0 = ha + (I, - 13) @3 @ (3) 
0 = 13 03 


b3 


= 


—e 


te 


Axisymmetnic Rigid Body 
Figure (1) 


The angular velocity vector, angular momentum vector, and the kinetic energy may be 
expressed respectively as 


N..B 


®O = @)b, + @2b2+ 3 b3 (4) 
h = 1]; @; b} +12 @2 bo +13 @3 D3 (5) 
E= Ln wo? + I, ws + 13 w?} (6) 


A simplification of the notation can be made. Let/; =/2=/, and 4; =/,. From the third 
line of Equation (3) it can be seen that the angular velocity about the spin axis b3 is 
constant, therefore @3 = @; . The transverse angular velocity components interchange 


between the b, and the b2 axes but the magnitude is constant, so one can let 


2 


®, = @, b, + @2 b2. Therefore, 


N.B -. 
® = @,+0, b3 (7) 
h = 1,0; +15 @s b3 (8) 
he = leo? +12 02 (9) 


eS Lr, w,? +], o 2) (10) 


Because the motion is torque-free, |h| is constant and h is fixed in space. Because there are 
no energy sources or sinks, E is also constant. Finally, the above conditions result in 
ie being a constant. 

An expression for the nutation angle may now be developed. The orientation of 
the body axes with re: -t to an inertial reference frame is desired to provide a measure of 
the body's dynamic behavior. Because the angular momentum vector is fixed in space, the 
nutation angle is defined as the angle between the body-fixed axis about which spin is 


desired and the angular momentum vector, and can be expressed in one of the following 











forms 
6 = cos"! ae = cos"! (13) = cos7! a = cos7! a ‘ (11) 
6 = sin7! [yore ete = sin7! [fp SB Bea P| = sin7! ca (12) 
_; fh, bi + he bd ee + (hed 
x 1 {21 Dy + 22 = 1 {41 @1 D1 + 12 @2 - 1 1 
a ian i ba nee Se e,) OO 





In these equations the first and second terms correspond to the general case with spin about 
the b3 axis, and the third term uses the simplified notation to describe an axisymmetric 
body. In the special case of spin about only one principal axis of an axisymmetric body, 
the angular momentum vector and the angular velocity vector will lie on the spin axis. With 
spin components on two or more principal axes, the three vectors will not be coincident, 


although they still will lie on the same plane. 


2. Single Spin Satellite Stability 

A torque-free, axisymmetric rigid body with the body axes coinciding with the 
principle moments of inertia will be stable about the axis of either the maximum moment of 
inertia or the minimum moment of inertia. To prove this, one begins with an arbitrary rigid 
body. The body is given the inital condition of steady angular velocity, @p, about a 
principal axis, and is then perturbed slightly. It is assumed that the angular velocities about 
the other axes are small, and are approximately the same order of magnitude as the 
perturbation (o = W2 = e). The system will be considered stable if the perturbation does 


not increase over time. Given an initial angular velocity with a perturbation, 


N..B ' , 
W® = @ by + @2 b2 +(Wo +e) b3, and given arbitrary inertias /;, J), J3, Euler's 


equations of motion can be written as 
O = 1, @ + (13 - 12) 2 (a + €) 
0 = Ip @2 + (I; —13) @ (@ + €) (14) 
0 = 13 3 +(I,-h) @ a 

The equations are linearized by neglecting terms of magnitude €?. Rewriting the equations 


by eliminating the terms € @, € @, and @, @ results in 


. I,-] 
a, = 2-13) a wo (15) 
-  _ (13 -h) 

= --——- @, @ 
W2 le 1 WO (16) 
Ww = Ot+ée=eE=0 (17) 


From Equation (17), one can conclude that €= constant. By differentiating Equation (15), 


and using Equation (16) to eliminate w, one gets 


@) Hi U3 ~ h) (3h) 2 @, = 0 (18) 
I, 12 


Similarly, from differentiating Equation (16), and using Equation (15) to eliminate @), one 


arrives at 


i + V3 — 12) (13h) «2 a = 0 (19) 


Ii 12 
These equations are identified as second order, linear, ordinary differenual equations with 
constant coefficients. The general solution for these differential equations are 


(20) 


K, elt + K>2 e-!yt 


K3 e'Y' + Kg e-'7! (21) 


@| 


@2 


I3 -—12) (13 -1 
= , (a= b)b eed (22) 
I 12 


If y is imaginary, @, and @» will increase without bound over time, and the motion will be 


where 


unstable. Stability is achieved if y is real. The first case occurs when the maximum 
moment of inertia is about the spin (b3) axis; then /3 > J; , 1/3 > Jy, and (1/3 — 12) (13 — 1,) > 0. 
In this case the inertia ratio /, / J,, defined as the inertia about the spin axis over the inertia 
about a transverse axis, is greater than one. The second case occurs when the minimum 
moment of inertia is about the spin (b3) axis; then /3 < J; , 13 < Jy, and (J3 — 22) (/3 —1,;) > 0 
as well. Here the inertia ratio is less than one. In both of the above cases 71s real and the 
motion is stable. However, if /3 is the intermediate moment of inertia about the spin (b3) 
axis, then (/3 — /2)(/3 —1,) <0, 7 is imaginary, and the motion is unstable. 

The previous model cannot be applied to a satellite since the assumption of a rigid 
body cannot be extended to the spacecraft. Structural elasticity, liquid propellant slosh, 
etc., Cause energy dissipation in an actual spacecraft. This spacecraft can be generalized by 
a quasi-rigid body with an unspecified energy damper mechanism. A priori, one can 
conclude that energy in the above system will dissipate until the minimum energy state is 
reached. The kinetic energy of the system with spin about the principal axis with maximum 


moment of inertia and spin about the principal axis with minimum moment of inertia can be 


written respectively as 


oigeg \ Imax about b3 
Pigg oe 
5. (23) 
E = 22 I min about b3 
2 I in 


The angular momentum 1s constant in the torque-free case. Thus the minimum kinetic 
energy state occurs when the rotation is about the axis of the maximum moment of inertia. 
Therefore, a quasi-rigid body is stable only when it is spinning about its major axis, with a 
corresponding inertia ratio that is greater than one. 

A relationship can be established between the time rate of change of the nutation 
angle and the energy dissipation of a quasi-rigid, axisymmetric body. One must assume 
that the angular momentum and moments of inertia of the quasi-ngid body do not change 
appreciably from a comparable rigid body. For the generalized model, with arbitrary 


inertias /; = /2, and /3, Equations (5) and (6) are substituted into Equation (12) to obtain 
aor? 
in gee A (2 1; E-h?} (24) 
I3-1, h2 


The time rate of change of the nutation angle is determined by taking the derivative of the 


above equation 


ee eee ee ae (25) 
sin (26) (3 - 11) h? 

where the only rate of change of energy E is attributed to the damping mechanism and is 

written as Ep to1ai to emphasize this point. Because Ep to:ai is negative, the nutation angle 

will decrease only if /3 is greater than /,;. This reaffirms the previous conclusion that a 


quasi-rigid body is spin stabilized only about the axis of maximum moment of inertia. The 


foregoing development is referred to as the energy-sink method. 


B. DUAL SPIN SATELLITES 

The logical progression from the single spin satellite was to incorporate a de-spun 
platform. This permitted the replacement of the low gain omnidirectional antenna with 
directional antennas for communication satellites, and a more capable spacecraft for 
scientific observation. A simple control system about the b3 axis would maintain the 
platform rotating at a constant relative rate with respect to the rotor ( and would usually 
have the platform rotate at the earth's rotational rate). Initially, the platforms were 
sufficiently small, and the overall dimensions of the satellite were such that the inertia ratio 
would be greater than one. For this type of satellite, the previously developed theory 
proved adequate. However, as satellites continued to grow in size, the launch vehicle 
shroud diameter became a constraint. In order to provide the size spacecraft needed to 
satisfy mission requirements and still fit within the shroud, a spacecraft with an inertia ratio 
of less than one (Is totat/ J: total < 1) would need to be built. From the previously 
developed theory, a spacecraft with an inertia ratio of less than one was believed to be 
inherently unstable. It was not until the development of the energy-sink theory for dual- 
spin, quasi-rigid, axisymmetric spacecraft containing a driven rotor, developed 
Simultaneously by V. D. Landon (unpublished work) and A. J. Ionllo [Ref. 1], that a 
spacecraft with an inertia ratio less than one was considered feasible. Several mgorous 
stability analyses using the equations of motion for specific dual spinners have been 
performed by P. W. Likins [Ref. 2], D. L. Mingon [Ref. 3], and others, to validate the 
energy-sink theory. The difficulty of a ngorous analysis is in accurately modeling all the 
forms of energy dissipation. A more general and practical approach was required to 
determine stability, and the energy-sink theory proved suitable. The development of this 


theory is as follows. 


The spacecraft, shown in Figure (2), is assumed to fulfill the following conditions: 


¢ Both the rotor and the platform are axisymmetric 
¢ Both the rotor and the platform are quasi-rigid 


¢ The damping mechanisms do not significantly alter the energy value, although 
the mechanisms will contribute to the energy rate 


¢ No external torques are applied 
¢ The only relative motion is spin about the b3 axis 
¢ The motor, driven by the control system, inputs just enough energy to exactly 


offset the shaft frictional losses, maintaining a constant relative angular 
velocity between the rotor and the platform 





Dual-Spin Quasi-Rigid Axisymmetnic Spacecraft 


Figure (2) 








The magnitude of the angular momentum and the kinetic energy of the dual-spin system can 


be expressed respectively by the following equations 
r_. r\2 ,\2 
h? = led w,? + (1, O,+Is Ws ) = Loe) ow,” +. (7, O+ 1s total Ds (26) 


rd ,2 pe - 
z= ij. I, feigi e+ Ioz+t [,@Ms = a l, total Oy suk Ts totai@s +1 1,07 +1,0, o (27) 
4 2 ys 2 Z 2 


where oO is the relative rotation rate of the rotor with respect to the platform. Although the 
equations actually represent a rigid body system, they are also applicable to the quasi-rigid 
system because of the above assumptions. If the damping mechanisms do make significant 
contributions to the energy of the system, then Equations (26) and (27) do not hold, and 
the energy-sink criterion will not apply. Additionally, the potential energy of the system 
(for example the energy stored in the spring of a mass-spring-dashpot damper) will not 
make a significant contribution to the total energy of the system. Therefore, for the 
remainder of the thesis, it is assumed that the kinetic energy of the system is effectively the 
total energy of the system (E = Berit Because the system experiences no external torques, 
angular momentum is conserved. Because the motor contributes no energy to the system, 
the time rate of change of the kinetic energy E is represented by only the quasi-rigid energy 
dissipation Ep jo1a!, and is negative. Differentiating the above two equations with respect to 


time, one obtains 
0 = 17 d+ UI, @s +15 Os) (I, Os +1, Os) 2 
= Lf potagt Or Or + \Is Ws tls Ws J\ls Ost ls Ds (28) 
E = Ep total = 1; total O1 O, +1. Ws Os BF I, Ws Ws (29) 
Eliminating the common term @, QO; by combining the above equations 


+1,@,O,+1,0 W, O, (30) 
| t total 


a Ep total = 


One may now define the inertial nutation frequency, Ao, the rotor nutation frequency, 


A, and the platform nutation frequency, A’, as follows 


_ 1, @s +15’ Ws 





DS (31) 
I; total 
A = Ap— a, = Ys titoral) Os + Is! ws G2)— 
I; total 
Die arn = I, ow, + (I, —I; ail Ds (33) 
Ts total 
The nutation angle for a dual-spin system is defined as 
h 3 b3 h3 + h3’ I; Os +15) Os” 
6 = cos! mores = eee! | | = cos! ae (34) 
[h| Ih Ih 
By imposing the condition 
ho > 0 (35) 
the analysis of nutational motion is restricted to the following region without any loss of 
generality 
0<a<F (36) 


Incorporating the nutation frequency terms, the equation for energy dissipation is written as 
— Ep total = Ep+ Ep’ = —1,A @s- 1,’ Xr’ OW,’ (37) 
Recalling the assumption that the rotor and platform are uncoupled about the b3 axis, we 


may incorporate Equation (37) into the reaction torques which tend to change the angular 





rates 
I, O, > _&b 
(38) 
Ep’ 
i Ws pee D 
y 
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Combining Equation (30), (37), and (38), one arrives at the transverse rate equation 


En Ep’ 

I; total ©: @; = Ag [f2 + ‘ (39) 
A WR 

Differentiating Equation (12) and substituting Equation (39) into it, the time rate of change 


of nutation angle as a function of energy dissipation rates is 





9 - 2 tot 2a ED E (40) 
sin (2 a) can Wo, 


The energy-sink equation for a single spin-stabilized body is obtained by letting A and A’ 
become Ap and Ep + Ep’ become Ep joiqi for a single body. The definition of stability 
requires the nutation angle @ to remain constant or decrease as a function of time, so that 9 
is zero or negative. Because the factors outside the parenthesis on the night hand side of 
Equation (40) are positive, the stability criterion for a dual-spin, quasi-rigid, axisymmetric 


system becomes 


(41) 





One of the following cases will guarantee stability 


1) A>O and A’>0 














2) A>O, A’<0, and |-2-|<|2 
4a hi 
3) A<0, A’>0, and |-2|>| £2 
ae 














A specific example would be the model of a typical communications satellite, a prolate 


_ dual-spinner possessing an inertia ratio of less than one. In general, the rotor nutation 
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frequency, expressed as 


= I, Ms +15’ Ws’ a Us —T total) Ms +15’ Ws’ (42) 
I; total F I; total 


would be negative because (/; — J; sora) is negative and @, >> @, if @,’ is rotating at the 
earth's rotation rate. Thus, energy damping in the rotor, Ep, would be destabilizing while 
energy damping in the platform, Ep’, would be stabilizing. It is from this result that 
satellites will have a damping mechanism placed on the platform to improve nutational 


stability. Such a damper is called a nutation damper. 


192 


Hf. PROBLEM DEFINITION 


A. OVERVIEW 

The existing energy-sink theory relies on several assumptions, perhaps the most 
important relating to the driven rotor. As previously stated, it has been assumed that the 
motor, driven by the control system, inputs just enough energy to exactly offset the shaft 
frictional losses, maintaining a constant relative angular velocity between the rotor and the 
platform. In actual systems, contrary to this assumption, the motor may add or remove 
energy from the system, depending on the dynamics of the spacecraft. Consequently, a 
revised energy sink stability theory, properly accounting for the energy contribution of the 
motor, is derived. The revised theory, based on the concept of core energy, will remain 
consistent with the existing energy-sink stability criterion. Continuing, an equation for 
nutation angle over time, as a function of core energy, is developed. Given a postulated 
energy dissipation function modeling the nutation dampers, structural elasticity, fuel slosh, 
etc., one can accurately predict the nutation angle behavior. Numerical simulation of D. L. 
Mingori's dual-spin, quasi-rigid, axisymmetric system containing a driven rotor [Ref. 3] is 
used to validate the revised energy-sink stability theory. The predicted nutation angle, 
based on this revised energy-sink theory, and the postulated energy dissipation function, is 
compared to the exact nutation angle of the Mingori system. By using a suitable postulated 
energy dissipation function, one can achieve excellent agreement between the predicted and 
the exact nutation angle. I. Michael Ross [Ref. 4] performed this analysis on a dual spin 
system with a damper on the platform only. The remainder of this thesis will use the same 


analysis, but on the Mingori system with dampers on both the rotor and the platform. 
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B. INCONSISTENCY OF THE ENERGY SINK THEORY 

There exists a contradiction between the existing energy-sink theory and the nutation 
angle derived from it. This will provide the motivation for developing a revised energy- 
sink theory and an alternative equation for nutation angle over time. From before, the 


existing energy-sink stability criterion can be expressed as 








(41) 
and the nutation angle for the dual-spin system was defined as 
h- b3 h3 + hs’ I, @, +1, Ws 
= ot ee ee -1 = —— 4 
@ = cos | hl | cos | hl | Cos | i (34) 


As previously stated, for a prolate dual-spinner (Js sorat / Jt totat < 1), energy dissipation in 
the platform is stabilizing and energy dissipation in the rotor is destabilizing. The angular 
velocity of the platform about the spin axis, @,, can be expressed in terms of the angular 
momentum and the kinetic energy of the system. Combining Equation (26) and Equation 


(27), one arrives at 


wo,’ = - I, 0 + I; o ———E— O* (It totat — 1s) (42) 
Is total s total Is total (7 t total — 45 total) 

Substituting this expression into Equation (34) results in 

@ = cos-! (+ V0) (43) 
where Q is represented as 
7) = (2B Sivas) I, total + [30° | ie total (44) 
2 (1 = Is att} Ts total — 1s total 
It total 


Initial conditions at t= 0 will determine the correct sign, with continuity considerations 


maintaining the sign for all of r>0. Additionally, no external torques are applied to this 
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system, resulting in |h| being constant. Differentiation of Equation (43) results in 


@ = ee 8 eet (4G) = ,1 1 b£p ee total (45) 
1-(+ 40)" dt sind VO h2 (1 _/s a 
t total 


From the definition of nutation angle, Equation (34), and the condition imposed on it, 
Equation (36), the positive sign must be chosen in Equations (43) and (45). An important 
observation is made at this tme. Choosing the positive sign will result in a positive rate of 
change in the nutation angle, indicating an unstable condition. The relative rotor spin rate is 
an independent variable, and is arbitrarily selected here as a constant value over time. 
Therefore, energy dissipation 1n a prolate dual-spin spacecraft will make the nutation angle 
increase, regardless of whether the dissipation is in the platform or in the rotor. This is not 
consistent with the stability cnterion of Equation (41). Thus, the existing energy-sink 


stability criterion contradicts itself. 


C. CORE ENERGY AND ENERGY DISSIPATION 
The existing energy-sink stability criterion does not properly account for any energy 
that may be provided by, or absorbed by, the motor. To accurately represent the system, 
the total rate of change of energy must be written as 
E = Ep toa+ W (46) 
where W is the rate of work due to the motor, and may be either positive or negative and 
the rate of change of energy due to dissipative elements can occur in either the platform or 
the rotor. Recall that the kinetic energy of the dual-spin system was expressed in Equation 
(27). If the work due to the motor torque as a function of time is written in analytical form, 
then the time rate of change of the energy of the system due to all dissipative elements, 


Ep total, can then be expressed solely in terms of the quasi-rigid parameters of the dual-spin 


system. With this expression, the condition that Ep joa $ 0 will result in the required 
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stability criterion. The difficulty arises in that to get the work due to the motor torque W 
(or W), one needs to know the exact dynamics of the dissipative mechanism. In the 
development of the modified energy-sink stability, Equation (46) will be used to determine 
the expression for Ep jorqi, and ultimately derive the revised stability criterion and nutation 
angle equation. 

Additionally, the existing energy-sink theory can be shown to be incorrect for both 
the case of total energy decreasing and for the case of total energy increasing. For 
example, if the rotor is mgid and total energy decreases, Equation (41) predicts that the 
system will be stable, but Equation (45) predicts that it will be unstable. Allowing the total 
energy to increase would reverse the conditions, but still show a contradiction between 
Equation (41) and Equation (45). 

A modification of the energy-sink stability theory and the associated expression for 
nutation angle is now presented. The development of the theory is from I. Michael Ross’ 
unpublished notes. The basis of the new theory is centered on the core energy of the 


system. As defined by Hubert [Ref. 5] 


Core energy is the total energy of the spacecraft minus the portion of the rotor 
energy that is due to the relative rotation between the rotor and the platform. It 
1S assumed that the mass, inertia, and motion of the damping device are 
sufficiently small that its energy is negligible relative to the spacecraft core 
energy. The damper will be treated as an undefined ‘energy sink’ for the 
purposes of the energy sink analysis. 


From the above statement one can define the core body as that body whose inertial 
dynamics are selected for analysis. Hubert defines the platform as the core body. The core 
energy is simply the rotational kinetic energy of a fictitious rigid body that possesses the 


inertia properties of the entire dual-spinner but moves in inertial space exactly like the core 
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body. For a dual-spin spacecraft the platform core energy is defined as 
, 72 72 
EC = Ln total of +h total op +h; total 03 = Li, total o, +41, total Ds (47) 


The center expression is for an arbitrary dual-spin spacecraft with the spin axis about the b3 
axis, and the right expression is for a dual-spin, axisymmetric spacecraft with the 


simplified notation. Extending this concept to the rotor, the rotor core energy is defined as 
Ec = a total ©, + bo total OF + ae total 037 = Lh, total ©," + Li, total @s* (48) 


Necessary to the development of the modified theory is what will be termed the 
Separation Axiom. This is when analysis is first performed with the rotor considered ngid 
and the platform quasi-rigid. Euler's equations are written for the rotor, and through 
manipulation, an equation is derived relating the torque on the quasi-ngid platform solely in 
terms of platform vanables. Then the platform is considered rigid and the rotor quasi-rigid. 
An equation is derived relating the torque on the quasi-rigid rotor solely in terms of rotor 
variables. These two separate equations are then combined and applied to a system in 
which the rotor and the platform may both be quasi-rigid. 

The case of a rigid rotor with a quasi-rigid platform is first analyzed. Because the 
rotor is rigid, the torque applied by the motor to the rotor is the net torque on the rotor and 


is determined from Euler's moment equations. For the case of the axisymmetric rotor, 
Tr = Trim = Is Ws (49) 
One can observe that Tpyy=—Triy , but Tp#Tp since the damping mechanism 
contributes additional torques to the quasi-rigid platform. The rate of work needed by the 
motor torque to maintain constant relative motion between the rotor and the platform is 
W = Tro = 1,00, = 1, 0(@ +6) (50) 


By substituting the platform core energy, Equation (47), into the kinetic energy expression, 
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Equation (27), kinetic energy may be expressed as 
E=Ec' +31, o* +1, 0, 0 (51) 
The above equation is differentiated to arrive at the time rate of change of the kinetic energy 
of the dual-spin system 
E=Ec +lootl, (oa, + Ws 6) = Ec’ +], o(a,’ fs 6) +1,Q,/ 6 (52) 
Substituting into this equation the rate of work done by the motor torque, Equation (50), 
one gets 
E = Ec’ +W+1,0,6 (53) 
Comparing this to Equation (46), it can be seen that 
Ep torat = Ec’ + Is @, G (54) 
Taking the derivative of Equation (47) to get the time rate of change of the platform core 
energy 
Ec’ = Ir totat 1 @ + Is total Ds’ Ws’ (55) 
and then substituting this into Equation (54), one arrives at the total energy dissipation of a 


rigid rotor, quasi-ngid platform system 


Ep totat = 11 total O: @ + Ts total Ws Ws +15 Ws O (56) 
=I, voral W, WD, + 15 We We + 1," Ws. wo, +1, W, 6 
AY AY 


Because the system has no external forces applied, it remains torque-free. Thus, Equation 


(28) can be used to eliminate the @, @, term and arrive at 


Ep total = (7 s Ds + I a «,’) 7 
t total 


(57) 
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Noting that the platform nutational frequency, Equation (33), can be written as 


fae a = Z +1; mil @, = (7; total — 1; total) WO, +150 (58) 
I, total I t total 
then Equation (57) can be written as 
Ep totat = (Is s+ 1s @s)(- A’) = Ag i{- 2’) (59) 
Referring to Equation (34), the nutation angle can be written as 
anton ana Sad (60) 


Taking the derivative and comparing it to the rate of total energy dissipation, Equation (59), 


the following relationship can be established 


ar aman = = Ep total (61) 
h A [hl 
Because the rotor is rigid, all energy dissipation will occur in the platform, such that 
ED = |h| e'sin 6 (62) 





Ad 


Equations (61) and (62) can be rewritten by including the motor torque, Equation (49), as 








; _ En’ 
_|h| O@sin @ = Tr+I1, @, = —2 (63) 
Xn’ 
Because Tp)y = — Tp /y (action - reaction pair), the final relationship is wnitten as 
E ( Y Se 3 
Tp = ——- +1, @s (64) 


a 


This equation describes the motor torque on the quasi-rigid platform as a function of 


platform variables only. It can be seen at this point that the classical analysis can be 
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achieved by assuming no torque is applied by the motor, resulting in 7 pjy = 0 and 
ee ree I,’ W,” (65) 


This analysis can now be performed on the system containing a quasi-rigid rotor and 
a rigid platform. It is observed that the selection of the body to be the rotor and the body to 
be the platform is completely arbitrary, and no physical distinction exists between the two 
bodies. Therefore, by analogy, the equation describing the torque on the quasi-ngid rotor 


as a function of rotor variables must be 
E 
TRIM = - +1, @, (66) 


Now let Tp;y and Tp y represent the motor torques on the platform and on the rotor 
respectively for the system containing both a quasi-ngid rotor and a quasi-rigid platform. 


Then the separation axiom would require the following two conditions 


. Ep ee 
T py = Tpjm = =r tl QO; (67) 


: E : 
line = TR = — Os (68) 


Once again, since the system is an action - reaction pair 


and then 
2 E c : ° 
2p =D - (1, &, +1,’ @,') = [hl @sin 0 (70) 


AK 
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The stability criterion dictates that 9 remain zero or be negative, thus 
Ep + za < 0 :) 
A WwW 
This is seen to be the existing energy sink criterion. Therefore, despite the presence of 
motor torque, the existing stability criterion still applies. Referring to Equations (67) and 
(68), the second term in the equations would represent the torque applied by the damping 
mechanism to the platform and the rotor respectively. Equations (67) and (68) state that the 
motor torque minus the damper torque will equal the net torque of the platform and the 
rotor respectively. 
In determining the revised energy-sink stability theory, the energy dissipation 
equation for the system with no energy contribution from the motor, Equation (37), must 


be rewritten to account for the motor torque. Thus 
Es Ep total t W = —I,A@s—15 0’ Ws” (71) 


Substituting in Equations (67) and (68), one arrives at 


E 


E ° 
Ep trait W = 12 — Ti 





Ep’ * 
+,’ —-Tpy 
E | 


i ean (72) 
= Ep ote Ep ai Trim +A gn, 
By assigning the motor torque values as 
Tu = Trim = -Tpm (73) 
then the above equation can be rewritten as 
E = EpiatW = Ept+Ep’ +Tu (a -a) (74) 


Because the total time rate of change of energy dissipation equals the rate of change of 


energy dissipation in the rotor plus that in the platform, the rate of work due to the motor 
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torque must be 


W = Ty(4’-a) = Tyo (75) 

Referring to the example from Chapter 1, the prolate dual-spin communications 
satellite, typically the rotor nutation frequency would be negative and the platform nutation 
frequency would be positive. For rotor spin-up Ty, > 0 and the motor torque would add 
energy to the system, and for rotor spin-down Ty < 0 and the motor torque would remove 
energy from the system. For the arbitrary system, the motor torque, Ty, and the sign of 
the term A’ — A (= 0), would determine whether the motor adds or removes energy from 
the system. 

The rates of change of the energies of the system may now be represented. Rewntng 


Equation (46) to determine the total energy dissipation of the system 
Ep total = E-W (76) 
The rate of work due to the motor torque can be expressed by combining Equations (75) 


and (68) (or Equation (67) ) to arrive at 


. , Ep’ 
W = Pol QW; 0 -- 2 o- laos ; (77) 





a 


The rate of change of the kinetic energy of the system as a function of platform core energy 


and system parameters was determined previously as 
E = Ec +l, o(a@,’ “ 6) +1,@,. 0 (52) 
Substituting the above two equations into Equation (76), the expression for total energy 


dissipation of the system may now be written as 


Ep twtat = Ec’ +1; ola,’ + 6) +1,@;, O- a. o-1,@,0 (78) 
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Because @,’ + O = @,, this equation can be further simplified 
: a ,- Ep 
Eptoaat = Ec +1s@s O- Q 0] (79) 
Comparing this equation, representing the system with a quasi-rigid rotor and a quasi-rigid 
platform, with Equation (54), representing the system with a ngid rotor and a quasi-rigid 


platform, reveals an additional term, — = o. The first term of Equation (79) represents the 


rate of change of core energy, with the platform as the core body. The second term will 
account for the change in energy associated with a change in the relative rotation rate of the 
rotor with respect to the platform. The final term accounts for the energy loss due to the 
quasi-rigidity of the rotor that is not represented in the platform core body expression. 
The case that will be analyzed is that which occurs when there is a constant relative 
rotation rate between the platform and the rotor. For the remainder of the analysis, let 
o= 0 (80) 


and the total rate of energy dissipation becomes 


: -, Eno 
Eptotat = Ec - (81) 


Further simplification can be achieved by noting that o = A’—A and from Equation (37) 


that Ep total = Ep + Ep’. Therefore 





Ep+Ep = Ec’ -Ep 2 = : (82) 
which reduces to 


Ep ,Ep _ Ec (83) 
Aes. 


A similar development can be performed using rotor core energy vice platform core energy. 
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The result is 


y A 


(84) 


There is an expected symmetry between Equations (83) and (84), due to the arbitrary 
assignment of one body of the system as the rotor, and the other body as the platform. To 


confirm the results of Equations (83) and (84), one must prove that when o is constant 


EY 
geal (85) 
A Kr’ 
Taking the derivative of the rotor core energy, Equation (48) 
Ec = I oral O; : + Is toral Os O, | | (86) 
= I; total @; O, + Is retare a o) (a,’ Be c) 
and similarly, for the platform core energy 
Ec’ = It total O: @; +15 total Ds’ Ws” (87) 
Substituting Equation (87) into Equation (86) and noting that o = 0, 
Ec = Ec’ +1s torat 0 Os. (88) 
Muluplying through by the platform nutation frequency 
Ec XV’ = Ec’ + Ts total x’ O Os” (89) 
and recalling that A’ =A +0, 
EcNv = Ec’ A + Ec’ O+1s worat XO Ws (90) 
which is the same expression as Equation (85) if it can be proven that 
Ec’ + 1s total nr’ Ws = 0 (91) 


Eliminating the transverse angular velocity of Equation (47) by substituting in Equation 
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(26), one arrives at 


h? - (7, Ws + 1,’ w,’F 


Ej a ky ae le 92 
c=5 a > !s total Os (92) 
Taking the derivative 
I 7 , , : + ~ - : 
Ec’ = si S Os I, = (1, Qs I, Qs + I, pas, Ww, w,’ Z total (93) 
t total t total 
and simplifying 


Os a,') as S total I t total Os QO,” 


: ‘ : a 
E (2 Os Ws t+IsIs Ws Os tIsIs Ws i. Is" (04) 
4 I; sotal 


Noting that /; torat = I; +,’ and recalling that a, = @,’ because o = 0, the above 


equation can be reduced to 
ee ah _ Os +1, = — 1s toral Ds ) @, (95) 
t total 
finally, invoking the definition of the platform nutation frequency, one arrives at 
Ec’ = —1s total i Ws (96) 
Therefore, Equation (91) holds and the revised stability criterion can be expressed as 
(97) 





A few remarks can be made concerning this stability criterion. The third expression is the 
existing energy sink stability criterion. This criterion must equal the stability criterion as a 
function of the rotor core energy which must equal the stability criterion as a function of the 
platform core energy. It can be seen that one no longer needs to know the energy 
dissipation rates in the platform and in the rotor to determine stability. By knowing or 


postulating the core energy over time, the stability of the system can be determined. 
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Continuing with the prolate dual-spinner example of Chapter 1, any one of the above 
expressions apply. For a stable system, the rotor core energy will be positive, and increase 
over time. Additionally, the rotor nutation frequency will be negative, resulting in a 
negative expression for the rotor core energy stability criterion. The platform core energy 
will be positive, but will decrease over time. The platform nutation frequency will be 
positive, resulting in a negative expression for the platform stability criterion. According to 
Equation (97), the platform and rotor stability critena must equal one another. From the 
numerical simulation of a dual-spin quasi-rigid axisymmetric system, the rotor and platform 


core energies as a function of time will be determined and graphed. 


D. NUTATIONAL MOTION 

The development of a modified expression for the nutation angle as a function of time 
may now be presented. The actual nutation angle of the system is defined as 8. The 
nutation angle as a function of platform core energy will be represented by 77, and the 
nutation angle as a function of rotor core energy will be represented by 7. The derivation is 
similar to the one previously given in this chapter, except that the total energy of the system 
has been replaced by the core energy of the system to eliminate the transverse angular 
velocity,@,. The derivation for the platform core energy will be shown. From before, the 


angular momentum of the dual-spin, quasi-rigid, axisymmetric system is 


a? , Z t 4 2 
h2 = 12,,.,0,2 + (I, 5 + 1g Os)” = Figeg O,2 + I O + Tp tosat Os) (26) 
Combining this with the platform core energy, and solving for the platform angular velocity 


about the spin axis 


,2 , , 
I; oral (Is total — 1; total) 03 — 21s torat ls 0 W3 +12 o*+2Ec I; worat— h? = 0 (98) 
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and 


Ey’ = I, 6 +V0’ It total ] (99) 
(7, total — I; total) I, total — I, total I; total 
where 
6” 5.2 Ex Is rua 1 — Js total) _ 2 {4s total Zs soual)( _2 is total [seal 12.02 (100) 
I t total I t total I t total I t total 


The initial conditions at t = 0 will determine the proper sign of Equation (99). As before, 


continuity will ensure this sign for all > (0. Substituting Equation (99) into Equation (11) 


= cos"! | Frpworme poser fi '0’) (101) 


t total ~ I S total 


I 03° +I1,0 

a _1 [is total Y3 s 
= COS Fssaat s+ Ls 0 
7 h 





The time rate of change of the nutation angle will determine the stability. Differentiating the 


above equation results in 





7 = Fc 1 | (102) 


A stable system would require that the lower sign be chosen for the radical, thus 


1! = cos-! i Tt total {el lo-vQ o’}} (103) 


t total — I, total 


When this sign selection is applied to Equation (52), one gets 


M3 (1; total — I, total) )+ I; Oa vO [pow tea (104) 


S$ total 
This can be rewritten as the well established dual-spin stability condition, as written by P. 


C. Hughes [Ref 6], 


(I; totat — Ts total) 03° +1, 0 < 0 (105) 


It is important to note that this analysis does not produce a contradiction to Equation (41). 
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In a similar manner, the modified nutation angle can be derived with respect to the rotor as 


n = cos-! (|-—ttsoat__] 1 LiJo- re) (106) 


I, total — I; total 


where 


9) |p Is wall - Ts total) _ p,2 As soial}(y (ets {Geeta 7262 (107) 
I ft total I t total I t total I t total 


Therefore, 1f one is given the core energy or the postulated core energy as a function of 
time, the nutational motion can then be determined. This leads to an extremely important 
conclusion. By determining a sufficiently accurate postulaton of the core energy of a dual- 
spin, quasi-rigid, axisymmetric spacecraft over time, one can predict the nutational 
behavior and the stability of that spacecraft. Numencal simulation of such a system will 


confirm this conclusion. 
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HI. DYNAMICAL EQUATIONS 


A. MINGORI DUAL-SPIN SYSTEM 

A specific model is required to validate the modified energy-sink stability theory. The 
development of the previous chapters was completely general in nature and applies to any 
dual-spin spacecraft with an arbitrary damping mechanism. D. L. Mingori's dual-spin 
system [Ref. 3] provided the needed model required to validate the proposed stability 
theory, Figure (3). Additionally, the complete non-linear equations of motion were 
presented by Mingori; however, the expressions for the important parameters in attitude 
dynamics, angular momentum and kinetic and total energy, were not presented in his paper 
and had to be derived before any analysis could be performed. 

The Mingori system is comprised of two symmetric ngid bodies, the lower which 
shall be defined as the rotor, and the upper body shall be defined as the platform. By 
convention, all terms refernng to the platform will be the same notation as that of the rotor, 
except that they will carry the prime mark. Both the rotor and the platform have coordinate 
axes fixed to the body and located at the respective centers of mass. The distance between 
the centers of mass 1s specified by L. The entire spacecraft has coordinate axes fixed to 
the spacecraft center of mass, denoted by the double prime coordinate axes, and rotating in 
the same manner as the rotor coordinate axes. The coordinate axes b3, b3’, and b3"are all 
collinear. The spacecraft center of mass will vary along the b3 axis as the point masses in 
the rotor and the platform oscillate. The only relative motion of the platform with respect to 
the rotor is angular rotation about the b3 axis. A motor driven by a control system 
maintains a constant relative rotation rate 0. The angle between a line parallel to b, and 


passing through the platform center of mass, and D1’ is represented by y = o7. 
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b, b', b", 





Mingori Dual-Spin Quasi-Rigid Axisymmetric System 
Figure (3) 
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Each body contains a mass-spring-dashpot mechanism. An actual spacecraft 
undergoes damping from various mechanisms. Undesired damping can occur due to 
structural elasticity and liquid propellant slosh. Nutation dampers are incorporated into 
spacecraft to improve the nutational motion. The dual-spin axisymmetric system with 
mass-spring-dashpot dampers cannot accurately model an actual spacecraft, but 1s used to 
illustrate and validate the theory. A description of the rotor is as follows. The mass- 
spring-dashpot mechanism can be modelled by a particle mass mm attached to a spring with 
constant k inside a tube filled with viscous fluid with damping coefficient c. The motion 
of the particle is constrained parallel to the b3 axis only. At rest, the particle mass lies 
along the rotor’s b; coordinate axis, at a distance a from the rotor's center of mass. Three 
balancing masses, m2, m3, and m4, each equal to the mass of the mass-spring-dashpot 
mechanism, are rigidly fixed a distance a on the bz, —b;, and —b> axes. These masses 
maintain the axisymmetry of the system about the b3 axis. Displacement of the particle 
mass 7m, is denoted by the variable z. A simplification in this paper of the Mingori dual- 
spinner system is the assumption of a massless spring-dashpot system. Thus the particle 
of the mass-spring-dashpot system, 7m, is the same mass as the corresponding three other 
balancing masses, m2, m3, and m4. The platform can be described in a similar manner, 
with all notation modified with the prime mark. 

The dual-spin system center of mass coordinate axes rotates at the same rate as the 
rotor coordinate axes, but is located along the b3 axis at a distance Z,,, from the rotor center 
of mass. This distance will vary as the particle masses are displaced. Relating the dual- 
spin system's coordinate axes to the rotor coordinate axes was arbitrary. Equivalent results 


would be achieved by selecting the platform coordinate axes instead. 
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B. ANGULAR MOMENTUM 

The angular momentum of the Mingori dual-spin, quasi-rigid, axisymmetric system is 
now derived from first principles. 

The angular momentum (moment of momentum) of a particle of mass m; located in 


body B is defined as 


N=i 


h,=r;xm; V (108) 
where h; is the angular momentum of the i th particle with respect to the system center of 
mass, Fr; specifies the location of the i th particle with respect to the system center of mass, 
and Ny = of iy is the absolute velocity of the i th particle with respect to the inerual 
reference frame N. Expressing the absolute velocity in the system reference frame 


N--i  Nz-ecm 


V = Vaan 1; (109) 

where B represents the reference coordinate axes of the system. For the following 

derivation, all displacements and velocities are referenced with respect to the rotor (bj) 
coordinate axes. The angular momentum vector is rewritten as 


N= N-B 
hoa V eas @® x ri} (110) 


Applying this equation to particle 1 with mass mm, and position 


r, =a b, +0 by + (z—zZem) b3 (111) 
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one alrives at 


(112) 
(2 - Zom) 0 —a(zZ- Zcm) Qw, 
0 a? + (z - Zem¥ 0 @, /)+ 
b 
wale ze.) 0 a? O; . 
hy =m b, 
0 —{(z - Zcm) 0 Vim x | 0 | 
(z - Zcm) 0 —a Vem y } +| —a(z—Zem) b; 
0 a 0 Geen 2 0 
For the ngid body, Equation (110) is applied to a differential volume at location 
raw = ub; + v bo + (w— Zcm) b3 (113) 
Integrating, one arrives at 
I, +M zé, 0 0 QO; 
0 [> +M Zen 0 @> + 
b, (114) 
0 0 [3 3 
hy = b2 
0 M Zcm 0 Vem x | 
—~Mizcm 0 0 Vien y b; 
0 0 0 a 
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Similar calculations are performed for the remaining seven point masses and the platform. 


For the rotor, one arrives at 


In+M 22,4 


27-22 2Zem 


DM + 4m = 
20 Aca 0 
27-22 2m 
0 2a7+422,4 
—maz 0 
+ 
0 
, | (M44 m) (zm) 


—-mZz 


[ 0 


(M +4 m) (zm) 
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—maz ae 
0 @2 
I3 +4 ma? ‘¢ 
\ 
b> (115) 
a 
0 V om x 


0 | V cms 


Similarly, for the platform 


(116) 
bhi’ + 4m’ = 
I)’ + M’ (l-zem)' + 
, [207 +4 (I- zen) + 0 ~ m’ a’ 2’ cos (a1) 
z+ 2 z’ (1 - zm) @ 1 
I,’ + M’ (1 —zem) + 
0 ] icy Apiamly - m’ a’ 2’ sin (ot) Or 
ii 2”? + 2 2’ (1—zem) 
@3 
- ma’ 2’ cos(o2) -m'a’z'sin(o2) In’ +4 ma” 
—m’ a’ ocos(ot) 2’ + m’ a’ sin(ot) 2’ he 
+ | —m’ a’ osin(o?) 2’- m’ a’ cos(o8) 2’ b2 
(13° +4 m’ a”)o \,./ 
0 _~(M +4 m’) (1 — zem) 0 Yyieee 
—m'2 
ip (M’+4 m’)(I- zm) 0 0 Vee 


+ m’ 2° 


| 0 0 0 || Ven: 


The angular momentum equations have terms corresponding to the velocity of the 
center of mass. Angular momentum, when taken about the center of mass of a system, by 
definition must be independent of the translation of the system's center of mass (with 
respect to an inertial reference frame). To verify that this occurs in the above equations, the 
equation for the position of the center of mass is substituted into the center of mass velocity 
terms, and then these equations are equated. If they have the same magnitude but opposite 
sign, they will cancel each other out and it will be proven that the system angular 


momentum is independent of the translation of the system center of mass. 
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One must prove 


Weal EF) 2 hy sel a 
Looking at the b; components, one finds 
—|(M’ + 4m’) (—1 + zom)—m’ 2] =[(M +4 m) zem — m 7] (118) 
Rewriting this, one obtains 
(M’+4 m')l +mz+m' 2’ =(M+4m+M'+4m’) 2om (119) 


The center of mass with reference to the b;, bz, b3 coordinate axes is 


2, Mi i+ mit z+(M’+4 alee a 
; m 
lem = eee eee ee (120) 
>) M; +m; M+4m+M’+4m’ 
i 
Substituting Equation (120) into Equation (119) results in 
M+4m')l +m’ 2’ 
(M’+4m’)l +mz+m'z =(M+4m+M'+4m’) mz+(M’+4 m)1 +m’ 2 (121) 
M+4m+M' +4 nm’ 
This reduces to 
(M’+4m')l +mz+m'27 =(M'+4m)I +mz+m' 2 (122) 


Therefore, the angular momentum along the b, axis is independent of the velocity of the 
center of mass. Verifying this along the other axes can be done in a similar manner. This 
proves that the angular momentum of the entire system about the system center of mass is 
independent of the translation of the system center of mass. The angular momentum of the 
entire system can be written as 


h=hysamt Nha’ 44m’ (123) 
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where 


DM + 4m = 
I1+M 22,.+ 
= * 0 —maz ” 
z7-22 2m 
In+M zen + (124) 
0 on 0 @2 
ii z7-22 2m 
—maz 0 I3 +4 ma? ,? 
bi 
0 “ 
+ |—maz bs 
0 
and 
(125) 
DM’ + 4m’ = 
Ih’ + M’(l-zem) + 
Riga + 4(l—-zm)y + 0 -m' a’ 2’ cos(o1t) 
27422’ (i- Zen) | 
Io’ + M’ (l-zem) + 
0 apelin -m’ a’ 2’ sin(o1) ? 
eee el 2) 
O3 
-m’' a’ 2’ cos(o2} —m’' a’ 2z'sin(ot) Ix’ + 4m’ a” 
b; 
b2 
—m' a’ ocos(a t) 2’ + m’ a’ sin(ot) 2’ |, | 
+ | -m'a’ osin(ot)2’- m’' a’ cos(at) 2’ ° 
(13+ 4 m' a*)o 
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C. ENERGY 
The kinetic energy, E, and the total energy, E;o,a), are now derived. From first 


principles the kinetic energy of a differential particle is written as 


Bam = ¥| yam P dm = 5{( Mum) - (My) | am (126) 


where the absolute velocity is defined from before 


N--d my: N-B 
Vv ” = v +Tfant @ Xam (109) 


The kinetic energy of the Mingori system is determined by integrating the differential 
kinetic energy over the rotor and over the platform, and summing the differential kinetic 


energy equation over the eight point masses to arrive at 


c= RL Um) (a Name SHG) (Nt 
oe 7 (127) 


xf Maaee (ogee awe 3 AAma)- (6D 


i 
Performing the steps on particle mm), the following is obtained 


(128) 


L(myem =. Nem) a (ey - ri) + (Mom) - Maem) ) + 
ey 57 
Bim bm Se 


—-cm P N=ecm si sos 
i=1 V5 +2("V agers r1))+2 (1 - Nari 
Substituting in the values, one arrives at 


e e 2 
Vin xt Vem y* Vin 2 (2 — zem) + (z— zemP 3 + 
a2 03 2 +(z-ZemP w?-2a (Z — 2cm) @) @3 + 
2 ; : 
my a? ws + 2(2 Vem 2—Zcm Vem z) + (129) 
Z M2 Vem x — Z7em M2 Vemx + 4 3 Vem y — 
Z@1 Vem y + Zem M1 Vemy — @ M2 Vem z 


E 


Np 


fp 
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This operation is then performed over the other three particles in the rotor. Performing 
these same steps on the body of the rotor, 


(130) 


(yom. Mem) a (Fy #1) + ( (Mem 1). (NGO x M1) ) + 


2 (MI. J+ 2 (MI. Wom xr) J+ 2 (Fr - GERM) 


Cry 
x 
It 


Substituting in values, and specifying the position vector of the differential particle as 


rau = u by + v bp + (W—Zem) bg (131) 
then 
(132) 
Vz,x+ Vimy + Ven 2+ Zam + Ww? WF —2 W Zem ws + 
Zam OF + v2 WF -2WV 2 Wz +2 V 2m W2 W3 + W2 W? + 
Si 22m 0% + u2 w2-2w Zem OF —2 UW @) @3 +2 U Zem W, W3 + 
oe 2 v2 w2 -2 vu @ 2 + U2 Ws + 2 (—Zem Vem 2) + dM 


W 02 Vem x — Zem @2 Vemx- V @3 Vemx— W 1 Vem y + 
Zem @1 Vem y + U @3 Vem y+ V @1 Vemz—U @2 Vemz 
2(—V cm @1 + UZcm W2) 


+ 


Ror 


Integrating over the limits of the symmetrical body, the kinetic energy of the rotor becomes 
Ven xt Vey + Vim: + tem + zen (or mE w3) + 


M I; w? + [5 w3 +13 ws _ (133) 
2 Zem O2 Vemx + 2 Zem @1 Vem y — 2 Zem Vemz 


These same steps are performed on the platform and on the remaining seven particles. 


These equations are then combined and simplified. 
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The total kinetic energy for the Mingori system 1s 


(134) 
> Miotat (Vern xt Vin y+ Van :} 3 
1h 0? +1, 03 +13 o3)+1 (1! of +h’ @ @5 + 13’ (@3 + oP) + 
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1(;7 +2 “(2 + ‘o)- 2 a’ z’ cos(Ot) @ @3 +2 a’ 2’ sin(or) a ws) + 
m(—az@)+ 
E = | m'(-a’7'cos(o1) a) +a’ 2’ sin(ot) @ - a’ 7 ocos(ot) @ — a’ 2’ Gsin(ot) @>} 
1 (ot+ of) (mzt+m zf-(mi+m 2} | 
"9 Mr 
M’+4m’\ 





1 (a? + w3)(M+ 4m)? vm 
1 (aj + w2)(M’ + 4m’) 12 rae | = 
ce an 
Mr - 


4 O 4 
? (w? + of) (M2 z’—-(M + = 
Mr 


mz(@? + @3} 





The total energy of the Mingor system can be determined by adding the kinetic 
energy of the system to the potential energy of the system. The only potential energy of the 
system is the energy stored in the springs of the mass-spring-dashpot damper. From 


Hooke's law, the system potential energy 1s written as 


U(z)=Lk2relez 2’ (135) 


The total energy of the Mingori system is then 


Eval = E+U = E+Lkz2+d v2” (136) 


D. MINGORIS EQUATIONS OF MOTION 


Mingoni's equations of motion [Ref 3] are written as follows 


A @, —(A — C) @) @3- Jy’ 60. +2Mr 600, + Mz C7 (a; + 2 3) + 


| ~2(£+ b)[z (a - @2 @3) +2 @] + . 
z|-2 €@\+22@ +2(@, - a w;)-a(a3 + 00) 
—2 (= h)| 2” (a, — @03) +7’ o;| + 
m’ \z'|-—2 ba +27 @) + 70, - a ,0)-a' cosw(ds + @,0)| + =0 
a’ siny {2 +z’ [(@3 + of — @3 | | 


A @, —(C —A) 0, — Jy 6, + 2M CC +Mz C7 (a+ 0, @3)+ 


(137) 
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m(1—p)z-m’ pz’ -ma(a,- a @3)- 
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In the above equations, the following relationships are used 


Mr = M+M'+4mp+4m) 





a M’+4my) 
= 
= mM mn 
P Mr p Mr 
I = (M+ Smo) bp a (+4 me) 
Mr Mr 
(138) 
C= (pz+p'z) y= or 


Az=il, +1) +2 m, a2 +2 mz a” +(M’'+4m,')(1- v2 


Zz 
C = 134+1;'+4m,a2+4m, a’ 


J3’ = 1,°+ 4 mp’ ag 


Chapter IV discusses how these equations are adapted for use in the numerical integration 


routine. 


E. NUTATION ANGLE 
The nutation angle of the dual-spin system can be determined by substituting 
Equations (124) and (125) into Equation (34) to arrive at 
(139) 
& maz—m'a’z’ cos(o7)} @) + 
(= m’ a’ 2’ sin(ot) @ + 
e@= cos Ps] = cosf s+ = cos-! (r, +4ma2+I13°+4m’ a’’} 3 


+ (74° +4’ a’) Oo 


Ih 
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F. CORE ENERGY 

The platform core energy and the rotor core energy were defined in Equations (47) 
and (48). The nutation angle as a function of core energy was then developed and is stated 
in Equations (101) and (106). These equation can be used to predict the nutation angle of 
the system as a function of time. Because the core energy as a function of time is not 
available for the prediction, a postulated core energy must be developed. The initial value 
of the actual core energy and the postulated core energy must match and 1s dictated by the 
initial conditions. Parameters of the postulated core energy must be selected to accurately 
model the actual core energy, to include the final energy state and the rate at which the it 
approaches the final energy state. Two models were considered. 

1. Exponential Core Energy Model 

An exponental representation of rotor and platform core energy as a function of 


time are expressed as follows 
Ec postulated (t) = (Eco — Ec final) e-") + Ec final (140) 


Ec postulated’ (t) = (Ec;’ — Ec fina’) e'-"*) + Ec final’ (141) 

The initial core energies, Ec,, Ec), are determined by the initial conditions of the system. 

The final core energies, Ec final, Ec final, aS well as the exponential factors 7, 7’, must be 

selected. The methodology for selecting these values is explained in Chapter V, Computer 
Analysis. 

2. Verhulst Logistic Core Energy Model 
The exponential model, as will be explained in Chapter V, has excellent 
agreement for stable conditions, but performs poorly for the unstable conditions. As an 


alternative model, the following first order differential model is selected for the rotor and 
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platform respectively. 


7 (Ec,) (Ec finat) 
MC postulated Y Ec, + (Ec finat — Ec,) e(-rt) sae 


(Ec. )(Ec finat’) 


ee (143) 
Ec, + (Ec final’ - Ec,’) e-") 


Ec postulated (1) = 


This type of equation was first introduced by P. F. Verhulst to model human and 
other populations [Ref 7]. It is often referred to as the Verhulst equation or the logistic 
equation. Although population dynamics appears to be unrelated to the stability of a dual- 
spin system, the behavior over time of this equation compared to the stable and unstable 
systems provides some insight. Figure (4) shows the Verhulst equation versus time for 


varying initial conditions. 


Verhulst Logistic Growth Model 





Figure (4) Verhulst Logistic Equation Versus Time 


It can be noted that if the initial value of the variable (in our case core energy) is 
greater than the equilibrium value, it will approach the equilibrium value in a manner similar 
to that of an exponential decay. For initial values that are less than but within one half of 
the equilibrium value, the variable will asymptotically approach the equilibrium value. If 
the initial value is less than one half of the initial value, the variable over time has a 
somewhat different shape as it approaches equilibrium. Initially the slope is very small, but 
increases to a maximum at about the one-half of the equilibrium position. After this 
inflection point the variable approaches equilibrium asymptotically. The value of No equals 
K is an asymptotically stable equilibrium. A value of No equals zero is an unstable 
equilibrium. If the value of No is slightly greater than zero, then N(t) will, as t 9, 
achieve the stable equilibrium of K. These types of curves will be utlized to describe both 
the stable and unstable dual-spin system. The initial core energies, Ec,, Ec,, would be 
determined by the initial conditions of the system. The final core energies, Ec final, 
Ec final » must be known or a best estimate used. The exponential factors, r, 7’, must be 
determined experimentally, and are a function of the system parameters and initial 


conditions. Chapter V provides additional details on the selection process. 
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IV. NUMERICAL SIMULATION 


A. NUMERICAL SIMULATION EQUATIONS OF MOTION 
The equations for the dynamical quantities needed for analysis of the dual-spin 
system were previously developed. Manipulation was required to make these equations 
suitable for numerical analysis, and is described below. The computer code is included as 
Appendix B. 
1. Mingori’s Equations of Motion 
The equations of motion for the dual-spin system are listed as Equation set (137) 
and (138). The Runge-Kutta numerical integration routine necessitated that the equations 
of motion be a set of first order differential equations. Mathematical manipulation was 
required to put them in this form. Variables representing groups of terms (A;, B;, Cj, etc) 
were introduced to simplify the manipulations of the equations and provide a suitable 
format for incorporating the equations into the computer program. Mingori's equations can 
be expressed in terms of these variables and the five tme-dependent variables of motion 
A26 @ + Az7 @3 + Arg 2’ + Arg = 0 
B25 @2 + B23 @3 + By3 Zz + Bz 7 + By =0 
C19 @ + Cs @ + C 3+ C11 = (144) 
D,z+Dz2z' + D3 a@+Ds,=0 
E, 2+ E27 +E3@,+ Es a) + E\9=0 
where the variables A;, B;, C;, Dj, Ej, F;, and Z; are defined in the Notation section. 


Clearly, these equations are highly coupled. Through a series of manipulations, including 


substitution and combining sets of equations to eliminate common variables, one can arrive 
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at a series of first order differential equations suitable for numerical integration techniques 








dz y 
-° 
-F,-28 F, +28 
Dp | SCO 
Zot Zo +2 
a> L D3 °"'Ds3 
a epee Z, +L 
D3 , _ D3 
- D> : D> (145) 
D — ee 
D; *" D4 
Dg Dg 
_F,—¥8 Dg 
2 D3 : F *D, 
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a ee 
Dp =D 
Zit+— 8 Z3+=1 
: D3 Sais 
(Es-Bis) z+ (Co Are , Cs Bai), = , Cio Azg , C5 Bag Cur 
day _ 4, - \C Bs Ar, _C Bos} C Ang C Bos C 
dt » C10 A27 _ Cs B23 
CAz C Bors 
gO) _ 4, = A272 G, 424 5°_ Arg 
dt ie Are Ar6 
da _ , = —223 4,213; _ B21 3 _ B26 
a pc Bos B25 B25 


The equations remain coupled, so the sequence in which the equations are numerically 
integrated is important. Because the equations for z and the equation for z’ are functions of 
z and z’, z and z’ must be integrated first. Similarly, 7 and z’ must be integrated before @3, 
@; before @,, and @ before @. In the computer program, the seven variables that 
describe the motion are stored in a seven column matrix where @ = y[1][i], @2 = y[2] [j], 


@3 = y[3)[j), z= y[4] GJ, z= y(5)fi), 2° =y[6) GJ, and z’ = y{7][j). The index j identifies the 
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matrix location for each set of saved variables over time. All other needed quanntes can be 
calculated by using the seven time dependent variables describing the motion, and the 
parameters of the dual-spin system. These additional dynamical quantities are stored as one 
dimensional vectors in the computer. 

2. Angular Momentum 

Angular momentum was derived in Equations (123), (124), and (125). The 

angular momentum 1s a vector, but for verification of the conservation of angular 
momentum, only the magnitude is required. Therefore, the three vectorial components for 
both the rotor and the platform are computed individually as H1, H2, H3, H1p, H2p, H3p, 
and then the magnitude of the angular momentum, hfj], is determined and stored. 

3. Nutation Angle 


The nutation angle is determined using the previously determined relation 


h P ’ 4 r 
@ = cos-! ae = cos-! ~ a | = cos"! eae (34) 


Since it is a function of angular momentum, it can now be calculated and stored as 





theta [j]. 
4. Energy 
The total energy of the system, Equation (136), is determined by the sum of the 
kinetic energy, as derived in Equation (134), and the potennal energy of the particle masses 


of the mass-spring-dashpot damper system, Equation (135). Itis written as 
Eital = E+U = E+Lkst+l er” (136) 


In the computer program the kinetic energy, ke [j], 1s the sum of the constituents of 
Equation (134), identified as T1 through T13. The potenual energy 1s not explicitly stated 
in the program, but instead is included in the calculation of the total energy, etotal {j. 
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5. Core Energy and Postulated Core Energy 

The platform and rotor core energy of the dual-spin system is defined in 
Equations (47) and (48) respectively. All of the variables are available, such that in the 
computer program the energies are readily calculated as Ecp [j] and Ec [j]. The postulated 
nutation angle may now be computed. Equation (101) indicates a sign must be selected 
depending on whether the system is stable or unstable. The calculation of nutation angle 
over time in Equation (34) will indicate stability. The computer program incorporates 
conditional statements to assign the correct sign in the postulated nutation angle expression, 
Equations (101) and (106). In the derivation leading up to the postulated nutation angle, o 
has been defined as the relative rotation rate of the rotor with respect to the platform, and is 
a positive value. This provides a positive contribution to the angular momentum vector, 
thus providing stability to the system. In Mingori's equations of motion, the reference 
coordinate axes are fixed to the rotor, resulting in the relative rotation of the platform in the 
counter-clockwise, or negative, direction. The postulated nutation angle equations were 
developed using the former reference frame. To compensate for the difference in the 
reference frames, the postulated nutation angle equations in the computer program have o 
replaced by —o. 

The postulated core energy as a function of time is computed for the platform and 
the rotor, using Equations (140) and (141) for the exponential model and Equations (142) 
and (143) for the Verhulst model. A postulated core energy model that accurately models 
the actual core energy required several iterations to find the proper values of the exponential 
factor. 

Once the postulated core energy models are computed, the Q’ and Q values, as 


defined in Equations (100) and (107), are determined. The postulated nutation angles, 
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etah [j], and eta [j], are computed using Equations (101) and (106) with the proper sign 
previously determined. 

The revised stability criterion, Equation (97), requires the time rate of change of 
the core energy and the nutation frequency. Since both of these parameters may vary with 
time, an estimate is made to determine if the modified stability criteria correctly predicts 
stability or instability. The time rate of change of the core energy is approximated by taking 
the difference of the final and the initial core energy values, and dividing by the time of the 
simulation. This quantity is then divided by the initial nutation frequency. This stability 
quantity is computed for both the rotor and the core, but is indicated in Table (1), Summary 
of Analyzed Cases, simply as a negative or positive value. 


B. COMPUTER PROGRAM 

Essential to the verification of the stability criteria is the computer program that 
integrates the dual-spin system equations of motion, calculates other dynamical quantities 
of motion, and graphs the results. The system used was a Sun SPARC 2 workstation with 
the computer code written in C. Intermediate graphics results were created using an in- 
house computer graphics program. Final graphics output was performed by sending 
output data to the Deltagraph graphics package. The sequence of steps of the computer 
code are explained below. The computer code is included as Appendix B. 

1. Initialization 

The main computer program is compiled, along with the header file ‘rkk.h’ and the 

function ‘derivs' immediately preceding it. The function ‘derivs' contains Mingori's 
equations of motion rewritten into Equation (145), suitable for numerical integration. 
Included in the header file 'rkk.h' is the numerical integration routine ‘rk4,' the adaptive step 
size function 'rkqc,' and the driver for the numerical integration routine ‘odeint.' Also 


included are functions for creating and freeing the vectors and matrices used by the 
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computer to store the data, and a function for error messages. 
a. Variables 
The variables required for both the numerical integration routine and 
subsequent calculations are defined as global variables before the main program. Variables 
used only within a specific function or only in the main program are defined in their 
réspective functions. Symbolic constants are also defined before the main program using 
the #define statement 
The time dependent variables of motion of the Mingori system, 
@, M2, W3, Z, Z, 2’, Z, are stored in a seven column matrix. All other parameters to be 
graphed are stored in one dimensional vectors. Storage in the vectors and the matrix 
permits retrieval by the graphics subroutine for plotting the vanable over time. 
b. Input File 
The main computer program scans the input file for the necessary parameters. 
All parameters of the Mingori system, /;, [, 13, M, m, a, k, c, 1;', [2’, 13’, Mum, a’, k’, 
c’, and L, are controlled with the input file. Also, the initial conditions for the time 
dependent variables of the system, @), @, @3, Z, Z, 2’, 2’, are specified. The length of 
time of the simulation and a variable determining the desired accuracy are specified. 
Finally, the exponential factor and the final energy for the postulated core energy functions 
are read. 
2. Preliminary Calculations 
After the input values are read in, initial calculations are performed prior to the 
numerical integration routine. All of the defined terms used by D. L. Mingon [Ref 3] in his 
non-linear equations of motion are computed, as well as definitions required for the core 
energy calculations, /s, Is’, Is totals Lt [1's Tt torat, aNd hrigia. 


An option is available to specify critical damping in either the platform or the 
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rotor. By using the number 1000.0 in the input file for the damping coefficient, the main 
program will automatically compute the coefficient required for critical damping of the 
mass-spring-dashpot system and then use this value in all subsequent calculations. 

Also computed is the factor ‘dxsav', used in determining when to save data. The 
steps between evaluating the equations of motion can become small, particularly when high 
accuracy 1s desired. The interval required for graphics resolution is not as restrictive. 
Accordingly, the variables are saved only if the step is greater than the previously saved 
step by the factor ‘dxsav.' 

3. Numerical Integration 

Mingori's nonlinear equations of motion are numerically integrated by the fourth 
order Runge-Kutta method, with adaptive step size control. The computer code used is 
based on the Runge-Kutta method listed in Numerical Recipes in C, [Ref 8]. The adaptive 
Step size control permits larger integration steps during smooth, well behaved portions of 
the functions, and smaller steps during the more irregular sections of the functions. The 
integration routine accuracy can be controlled by a variable in the input file. 

Mingoni's non-linear equations of motion are contained in the function ‘derivs.' 
As described previously, the equations have been rewritten as a series of first order coupled 
differential equations suitable for numerical integration. 

The output of the time dependent variables of the system, @), @2, @3, 2, Z, 2’, Z’, 
is Stored in an array of seven columns, with the number of rows required a function of the 
specified time interval and accuracy. A one dimensional vector is also created, with the 
time stored for each step saved. This permits plotting the time dependent variables and 
other quantities versus the time. 

4. Calculation of System Parameters 
With the array of the time dependent variables of motion as a function of time, the 


other system quantities listed in Section A may now be calculated, with the values stored in 
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vectors suitable for graphing. 
5. Graphics Output 

The remainder of the program is the necessary code for the in-house graphics 
program. The graphics program is initialized, and the graphics window is opened. The 
following parameters are then plotted over time: @, @2, @3 — @3 initial, Z, Z, 2’, Z , 
h — initial, E — Ejnitiai, Erotat — Etotat initia Ec — Ec initial: Ec postulated — Ec postulated initial: 
E~ — E¢ initial» Ec postulated — Ec postulated initial, 9, 7, and 1. The window is then closed 
and the initial conditions and pertinent time dependent variables and other quantities of the 
simulation are then printed. This displays the graphical representation of the behavior of 
the dual-spin system, and provides the actual initial and final values of the variables and 
associated quantities. 

For final graphics output, the Deltagraph graphics program is utilized. Computer 
simulation data is sent as an output file. The data is then manipulated into a suitable graphic 
with proper scaling and axes limits to best represent the dynamics of the computer 
simulation. The graphs are contained in Appendix A. 

6. Computer Program Validation 

The two aspects of the computer code requiring validation were the numerical 
integration routine and Mingori’s equations of motion with its associated dynamical 
quantities for the dual-spin system. 

The validation of the numerical integration routine was performed by using the 
equations of motion for a simple torque-free axisymmetric body. Various initial conditions 
were read in, and the dynamical quantities were then plotted versus time. The plotted 
behavior was then compared to the actual response of the system. 

The validation of Mingori’s equations of motion and the associated dynamical 


quantities required a more comprehensive approach. The equations of motion required 
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validation for proper behavior of each time-dependent variable of motion. Then the kinetic 
energy, total energy, and angular momentum must be verified. Several cases were run 
where the platform mass and inertia would become infinitesimally small. Then cases were 
run where the rotor mass and inertia would become infinitesimally small. In each of these 
cases, subcases were run where the mass-spring-dashpot system would be large, 
dominating the dynamics, to the subcase where it became very small, so the system 
approximates a ngid body. These different scenarios would uncouple and isolate the 
various parts of the equations of motion to verify proper derivation of the equations. In 
general, the platform and the rotor were tested under conditions varying from those in a 
rigid body scenario to those in a lightly damped body. The system was then tested as a 
rigid dual-spin system, as a dual-spin system with rotor damping only, and as a dual spin 
system with platform damping only. The dual-spin system was then tested with both the 
rotor and the platform containing dampers. For each case, the tme dependent variables of 
motion were plotted and analyzed. In all cases, the angular momentum was compared with 
the initial angular momentum. Since all cases were torque free, the angular momentum 
must remain constant. In all cases explored, the angular momentum verified the 
correctness of the equations of motion. The code validation cases were not included in this 


thesis due to the large number of graphs and data that were required to establish validation. 
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V. ANALYSIS 


A. INTRODUCTION 
1. Objective 

The development of the revised energy-sink stability theory has been presented in 
Chapter II, and the equations of motion for a dual-spin system are contained in Chapter II. 
With the numerical integration code of Chapter IV, the revised energy sink stability theory 
can now be verified. The core energy of the system is plotted as a function of tme and 
compared to the total energy for agreement with theory. The stability criterion computed 
from the numerical simulation must then agree with Equation (97). Conservation of 
angular momentum is used in each case to verify the correctness of the equations, and to 
ensure the accuracy of the numerical integration routine. To approximate the core energies 
Over time, an exponential model and a model based on logistic growth are then explored. 
The logistic growth model uses an equation presented by Verhulst [Ref 7], and is referred 
to as the Verhulst model. Postulated models of core energy and the associated nutation 
angles are then compared with the actual core energies and nutation angle for agreement. 

2. Numerical Simulation Cases 

Four distinct cases needed to be addressed. With the inerua ratio of the dual-spin 
System greater than one, a stable case and an unstable case are analyzed (Cases (1) and 
(4)). With the inertia ratio less than one, a stable case and an unstable case are also 
analyzed (Cases (2) and (3)). For these four cases, an exponential model is used to 
postulate the core energy and the corresponding nutation angle. For the case of the inertia 
ratio greater than one and unstable, and the two cases of the inertia ratio less than one, 


stable and unstable, the Verhulst model is postulated (Cases (5), (6), and (7)). The case of 


= 


the inertia ratio greater than one and stable was not included; from the excellent agreement 
of Cases (5), (6), and (7), one can see that excellent agreement would also occur for this 
trivial case, making it unnecessary to include. All seven cases are summarized in Table (1). 
For each case, the dynamical quantities are plotted versus time until the system reaches a 


stable state. Each quantity is then plotted for the first one hundred seconds to show detail. 


, 


T's total : EC’ 
CASE !twtal STABILITY x? MODEL COMMENTS 


large damping 
in platform 
large damping 
in platform 
insufficient damping 
in platform 
insufficient damping 
in platform 
same conditions 
as Case 2 


1.039 stable Negative negative exponential 
0.907 stable negative negative exponental 
0.897 unstable positive positive exponential 
1.025 unstable positive positive exponential 


0.907 stable negative negative Verhulst 


Methulst same condi 


0.897 unstable positive positive aa 


same conditions 


1.025 unstable positive positive Verhulst as Case 4 





Table (1) Summary of Analyzed Cases 


In Apper “x A, tables (2) through (8) lis’ xe system parameters, core energy 
parameters, and : conditions. Immediately follc «ing these tables are the graphs of the 
important dynar:. quantities. A typical geosynchronous dual-spin satellite was selected 
for the numerical simulation. The platform and rotor masses and inertias are listed as part 
of the initial conditions, and are the same for all cases. The inertia ratio is made greater 
than or less than one by selecting L, the distance between the rotor and platform centers of 


mass, to be 0.3 or 1.0 meter respectively. As a default set of values, all of the mass- 
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spring-dashpot systems have m and m’ equal to 1.0 kg, & and k’ equal to 1.0 N/m, and c 
and c’ equal to 1.0 kg/sec. To establish the stable cases, additional energy dissipation is 
required in the platform. To achieve this, cases (1), (2), and (5) have m’ increased to 20.0 
kg and c’ increased to 10.0 kg/sec. In all seven cases, the initial conditions are the same. 
The platform rotates at a geosynchronous angular velocity, the rotor spins at the higher rate 
of 1.5 rad/sec, and a perturbation is introduced by an initial transverse angular velocity of 
0.1 rad/sec. The mass-spring-dashpot systems have no initial displacement or initial 
velocity. The core energy parameters are also listed in the tables. The initial core energies 
are determined by the system’s initial conditions. The final core energies were taken from 
the numerical simulation data. The exponential factors, 7 and 7’, were then determined 


through an iterative process to best fit the modeled core energy to the actual core energy. 


B. DISCUSSION 

The individual cases can now be analyzed. The graphs of the dynamical quantities 
are explored, and the data will affirm the revised stability theory. 

1. Angular Momentum 

For each case, the angular momentum is plotted versus time, Figures (5), (6), 

(16), (17), (27), (28), (38), (39), (49), (50), (60), (61), (71), and (72). Because the dual- 
spin system has no external forces, angular momentum must be conserved. For the stable 
cases, there is excellent agreement, with angular momentum varying by less than one one- 
hundredth of a percent over the length of the data run. This confirms the equations of 
motion and validates the accuracy of the numerical integration routine. For the unstable 
cases, the angular momentum percent difference increases to approximately four one- 
hundredths of a percent. This is attributed to the increased dynamics of the system as it 


establishes the spin about the transverse axis, introducing very small errors in the numerical 
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integration routine. The errors remain very small, and the equations of motion and the 
numerical integration routine remain accurate. 
2. Total Energy, Platform and Rotor Core Energy 

The graphs of energy are in Figures (7) through (10), (18) through (21), (29) 
through (32), (40) through (43), (51) through (54), (62) through (65), and (73) through 
(76). Several observations can be made. 

The total energy curve reaches equilibrium before both the platform and the rotor 
core energies do. For Case 2 and 5, the total energy appears to reach a maximum and then 
drop down to a final equilibrium value. Careful comparison of Figures (18) and (49), 
shows the same system with the same initial conditions, but with a slightly different curve. 
It can be deduced that the energy is not steady, but 1s still oscillating. The sampling 
frequency coincidentally saved values near the same magnitude of energy in the region 
from about 3000 seconds to 10000 seconds. 

The total energy for Cases (1), (3), (5), and (6) decreased. For Cases (2), (4), 
and (7), representing both stable and unstable systems, it increased. This can be explained 
easiest with the use of equations (26) and (27). For Cases (2), (4), and (7), the system is 
settling out about the axis with the minimum moment of inertia. Because there are no 
external forces, the angular momentum is constant. In the equation for angular momentum, 
Equation (26), the inertia terms are squared. But in the equation for energy, Equation (27), 
the inertia terms are not squared. Therefore, as the angular velocity transfers to the axis 
with the minimum moment of inertia, Equations (26) and (27) show that the energy will 
increase. These equations apply to a simple dual-spin system. Equations (123) and (136) 
for the Mingori dual-spin system would show the same result, provided the energy 
absorbed by the mass-spring-dashpot system and the energy associated with the motor is 
less than the energy increase associated with transferring the spin to the axis of minimum 


moment of inertia. This is the situation for Cases (2), (4), and (7). 
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The total energy value will always be between that of the platform core energy 
and the rotor core energy. The platform core energy will have smaller values as its 
equation assumes that the rotor is rotating at the same rate as the platform, thereby not 
accounting for the large amount of kinetic energy associated with the rotor. The rotor core 
energy will be higher than the total energy as it assumes that the slowly spinning platform 
is rotating at the same rate as the rotor, providing the system with additional kinetic energy. 

The hundred-seconds graphs of total energy versus time shows curves with 
several different behaviors. The time rate of change of total energy includes the energy 
dissipation rate of the rotor and platform mass-spring-dashpot systems, and the rate of 
work due to the motor torque maintaining the relative rotation rate. In the hundred-seconds 
graphs, aside from the general slope of the curve, there is no apparent correlation between 
the specific behavior of the total energy to that of the core energies. Any relationship that 
exists is masked by the motor and mass-spring-dashpot system's influences on total 
energy. 

3. Stability Criterion 

The revised stability criterion of equation (97) states that the time rate of change of 
the core energy over the respective nutation frequency must be less than or equal to zero. 
The sign of the stability criterion for each case must be determined. The time rate of change 
of the core energy was determined by dividing the final less the initial value of the core 
energy by the length of time of the case. The nutation frequencies were computed using 
initial conditions. The sign of the stability criterion is listed in Tables (1) through (8). The 
numerical value was not listed because the above assumptions used during the computation 
give it no merit. Cases with inertia ratios very close to one were intentionally selected to 
test the inertia ratio in this transition region. In all cases analyzed, the sign of the revised 


Stability criterion was consistent with the stability of the system. 
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4. Postulated Core Energy and Nutation Angle 
An explicit relationship for core energy as a function of time does not exist. If an 
equation describing core energy did exist, then by using Equations (101) and (106) the 
nutation angle as a function of core energy and time could be predicted for a dual-spin 
system. This could then be extended to an actual dual-spin satellite. Given a sufficient 
model of the satellite, the stability and the nutation angle as a function of time could be 
predicted. An objective of this thesis is to see if an equation for the core energy can be 
developed that would adequately describe the nutation angle as a function of time for both 
the stable and unstable conditions. The exponential model and the Verhulst logistic model 
were explored. 
a. Exponential Core Energy Model 
Observation of the core energy as a function of time for a stable system would 
lead one to conclude that it behaves in an exponential manner. Equations (140) and (141) 
are exponential representations of rotor and platform energy as a function of tme. The 
initial core energy was determined by the initial conditions of the dual-spin system. The 
final core energy was determined by running the numerical simulation and calculating it at 
the end of the simulation. If this was not available, a value could be estimated. Applying 
the principle of conservation of angular momentum and noting that the system will 
eventually spin about one of the primary axes, then the equations for angular momentum 
and total energy can be used to solve for the angular velocity about that axis. Substituting 
this into Equation (47) or (48), one would arrive at an estimated final core energy. 
Although it would not be the actual final core energy, as motor torque contribution and the 
mass-spring-dashpot system’s energy dissipation was not accounted for, it would be 
sufficiently close to satisfy the computational requirements. The exponential factor is 
dependent upon the system parameters and initial conditions. For the cases presented, 


different values were tried until good agreement was established with the core energy 
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curve. Additional research would be required to determine a suitable exponential factor for 
an actual satellite, taking into account the parameters and the initial conditions. Cases 1 
through 4 contain the numerical simulation data for the exponential model. The actual and 
postulated curves for core energy are contained in Figures (11) through (15), (22) through 
(24), (33) through (37), and (44) through (48). For cases 1 and 2, both stable, there is 
excellent agreement between the core energy and the postulated core energy. This in turn, 
results in excellent agreement between the actual nutation angle and the modeled nutation 
angle. The actual nutation angle is determined using the angular momentum quantties, as 
shown. =quation (34). The time dependent variable required to compute the nutation 
angle is the angular velocity about the spin axis. The modeled nutation angle is determined 
in Equations (101) and (106), and 1s a function of only one time-dependent variable, core 
energy. Herein lies the potential of the core energy theory. A sufficient model of core 
energy over time, as in Cases 1 and 2, will provide an excellent prediction of nutation 
angle, without requiring any knowledge of the specific angular velocities of the system as a 
function of time. 

Cases 3 and 4 illustrate the exponential model for an unstable dual-spin 
system. The exponential model for core energy, and its associated nutation angle, rapidly 
approach their final values. The actual core energy and nutation angle, however, behave 
quite differently. The initial conditions have the system near an unstable equilibrium. The 
system moves from the unstable to the stable equilibrium slowly at first. It then increases 
the rate at which it approaches stable equilibrium, passes through an inflection point, and 
then approaches equilibrium asymptotically. It is clear that the exponential model 
represents this behavior poorly. Verhulst’s logistic equation was then addressed to 


determine its adequacy in modeling the dual spin system. 
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b. Verhulst Logistic Core Energy Model 

The Verhulst logistic equation was introduced in Chapter I]. Figure (4) and 
the associated description explains the characteristics of the equation. The general shape of 
the curves in Figure (4) is very similar to the stable and unstable cases of the dual-spin 
system. Case 5 is the Verhulst model of the Case 2 stable system. Figures (55) through 
(59) show that the Verhulst logistic equation can achieve excellent agreement with core 
energy and nutation angle, just as the exponential model did. 

Cases 6 and 7 are the same as Cases 3 and 4, except the Verhulst logistic 
equation is used to model the core energies. Figures (66) through (70) and (77) through 
(81) Ulustrate the modeled and actual core energies. Although the Verhulst model for rotor 
core energy had excellent agreement, it performs poorly when modeling the unstable 
system. The rotor core energy begins at an initial value and then decreases to an 
equilibrium value. Referencing Figure (4), the Verhulst logistic equation will model the 
rotor core energy exponentially. 

To take advantage of the Verhulst logistic equation’s curve beginning near the 
the unstable equilibrium and its progression to the asymptotically stable equilibrium, the 
core energy must increase over time. The platform core energy behaves in this manner as 
the nutation angle goes from a small angle to ninety degrees. Figures (68) through (70) 
and (79) through (81) show the Verhulst modeled and the actual platform core energies, 
and the modeled and the actual nutation angles. The agreement between the actual and 
modeled energies and nutation angles was very good. The general shape of curve was 
consistent, with only a slight deviation in the center of the curve, and then a small deviation 
as the nutation angle approaches the equilibrium value of ninety degrees. With this 
agreement, it is established that a core energy model exists that can represent both stable 


and unstable cases. For each of the cases, core energy 1s plotted versus nutation the angle, 
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Figures (15), (26), (37), (48), (59), (70), and (81). The graphics package permitted only a 
semi-log plot of rotor core energy. The plot did not provide any insight on relationships 
for the dynamics of the dual-spin system. For both the stable and unstable cases, the log- 
log plot between platform core energy and the nutation angle is linear from ninety degrees 
until about two degrees. For nutation angles less than two degrees, the platform core 
energy asymptotically approaches the equilibrium value. Case 5 and 6 is the same system 
with the same initial conditions, with the exception of the increased damping system in the 
platform for Case 5. The initial platform core energies and nutation angles are very nearly 
the same for each case, as shown in Figures (57), (58), (59), (68), (69), and (70). One 
can conclude by this observation and the definition of nutation angle as a function of core 
energy, Equation (101), that there exists a continuous platform core energy versus nutation 
angle curve. The appearance would look similar to the curve created by splicing Figures 
(59) and (70) together. By varying one parameter, for example platform energy dissipation 
rate, the system would progress along this curve and achieve equilibrium with a zero 
degree nutation angle or achieve equilibrium with a ninety degree nutation angle. This is 
what was done with cases 5 and 6. Adding as the third dimension to the curve the time rate 
of change of core energy would then reveal the stability, the initial direction, and the rate at 


which the system will arrive at the equilibrium condition. 


C. FURTHER RESEARCH 

The revised stability criterion for a dual-spin, quasi-rigid, axisymmetric system was 
established. Numerical simulation was then used to verify the revised stability theory. 
Further research could be conducted in several areas. The specific contributions to the total 
energy could provide some insight. How much energy and in what manner does the motor 
torque contribute to the total energy for both the stable and the unstable cases? Also, by 


plotting the energy dissipation system’s contributions over time, one could determine its 
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effect on total energy, as well as comparing it to the core energy theory in Equation (97). 
Relationships for core energy as a function of time were explored. It was established 
that the Verhulst logistic equation could be applied to a stable dual-spin system. By 
changing its parameters, this equation also applied to the unstable dual-spin system with 
very good agreement. Further research could be directed to find one equation that would 
address both the stable and unstable cases of the dual-spin system without changing its 
parameters. The equation for logistic growth with a threshold [Ref. 7] as applied to the 
platform core energy shows promise in this regard. 
Ec’ r __ EC 
Ec unstable Ec final 
Given the initial conditions, the platform core energy will typically start at some 





jae +(i = ge (146) 








intermediate value and will either increase or decrease as the dual-spin system reaches 
equilibrium with a nutation angle of either zero degrees or ninety degrees. Equation (146) 
would be well suited to model the platform core energy. There exists a platform core 
energy value, Ec unstable, Where the system is at an unstable equilibrium. Equation (146) 
shows that at exactly this value, the rate of change of the platform core energy will be zero. 
For any value below the unstable equilibrium value, the core energy will approach the value 
of zero. For Case 5, the final core energy for the stable condition was 0.122 J. Although 
this final condition cannot be represented in Equation (146), it is sufficiently close that the 
equation may still be used to represent the platform core energy. For any value above the 
unstable equilibrium value, the core energy will approach the equilibrium value of Ec fina! » 
associated with the nutation angle at ninety degrees. Initial conditions would determine the 
initial core energy. The final core energy can be estimated as described earlier in this 
chapter. Finally, the exponential factor r’ may then be determined based on the system 
parameters and initial conditions. 


Additional analysis can be performed on the revised stability criterion. The time rate 
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of change of core energy over nutational frequency may be calculated and plotted versus 
time or versus other parameters. The behavior of this stability criterion and the magnitude 
of it for various conditions could provide some insight on postulating a relationship for 


core energy as a function of time. 
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VI. CONCLUSION 


The existing energy-sink stability criterion was introduced as 


(41) 





It was shown that an inconsistency in the development disproves the existing 
assumption that the motor energy input exactly balances shaft frictional losses. A revised 
energy-sink stability criterion was then developed based on Hubert’s definition of core 


energy and was presented as 


(97) 





This criterion compliments the existing theory. Numerical simulation was required to 
validate the theory. The Mingori dual-spin, quasi-rigid, axisymmetric system was selected 
for the numerical simulation. Several cases were analyzed to verify the revised energy-sink 
stability criterion. By correctly postulating the platform or the rotor core energy, the 
stability of the system could be determined. Specific knowledge of the energy dissipation 
rates for both the platform and the rotor are no longer required. 

An exponential model and the Verhulst logistic model for core energy, and their 
relationship to the nutation angle, were explored. The exponential model had excellent 
agreement with the stable cases, but was inadequate in representing the unstable cases. The 
Verhulst logistic model established that an explicit relationship for core energy could be 
developed. Nutation angle as a function of core energy for the dual-spin system could then 
be predicted. The excellent agreement of the postulated core energy and nutation angle with 
the actual core energy and nutation angle confirms the revised energy-sink stability 
criterion. Additional research is required to find an optimum equation for core energy as a 


function of time to represent both the stable and unstable cases. 
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APPENDIX A NUMERICAL SIMULATION DATA 


CASE 1: a 1, Stable, Exponential Model 


f 


System Parameters Ininal Conditions 


Platform Rotor Platform Rotor 


0.0 m 


1,’ = 1600 kgm? 1000 kg m2 
0.0 


= 1500 kgm? 1200 kgm2 
1000.0 kg 
20.0 kg 1.0 kg 
1.0 m = 1.0 m 
N = N 
1.0 = 1.0 = Core Energy Parameters 


(fl 5 rad 
@3 = 7.27 x 107 <6 





= rad 
@), = 0.10 oe 


kg Platform Rotor 


kg 
tag = 0 ie 


Eni = 13.402 J E¢ initial = "OS Vane 








fstotal = 1.039 | |Ecfina’ = 0.0723 = Evfinat = “SlOL7G 
t total E-’ E 
= = negative —" = negative 
nr’ A 
r = -.00134s-! r = -.00134 s- 





Table (2) Case 1 Parameters 
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Case 1: Is/It>1, Stable, Exponential Model 
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Figure (5) Total Energy and Percent Difference Angular Momentum Versus Time 
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Case 1: Is/It>1,Stable, Exponential Model 
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Figure (6) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 1: Is/It>1,Stable, Exponential Model 
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Figure (7) Total Energy and Rotor Core Energy Versus Time 





Case 1: Is/It>1, Stable, Exponential Model 


energy (J) 
core energy (J) 


Ee 
6 on 
6 8 


woote* 
"ec 
e 
e 


grvgeenee” 
Peco le 


=" Ree., 
@e. 





Figure (8) Total Energy and Rotor Core Energy First 100 Seconds 
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Case 1: Is/It>1, Stable, Exponential Model 
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Figure (10) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 1:Is/It>1, Stable, Exponential Model 
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Figure (12) Modeled and Actual Rotor Core Energy and Nutation Angle- First 100 Seconds 
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Case 1: Is/It>1, Stable, Exponential Model 
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Figure (13) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 





Case 1: Is/It>1, Stable, Exponential Model 
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| Figure (14) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 1: Is/It>1, Stable, Exponential Model 
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Figure (15) Core Energy Versus Nutation Angle 
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CASE 2: is <1 , Stable, Exponential Model 
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System Parameters Initial Cond:ic as 
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Table (3) Case 2 Parameters 
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Case 2: Is/It<1, Stable, Exponential Model 
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Figure (16) Total Energy and Percent Difference Angular Momentum Versus Time 





Case 2: Is/It<1, Stable, Exponential Model 
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Figure (17) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Figure (18) Total Energy and Rotor Core Energy Versus Time 





Case 2: Is/It<1, Stable, Exponential Model 
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Figure (19) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 2: Is/It<1,Stable, Exponential Model 
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Figure (20) Total Energy and Platform Core Energy Versus Time 





Case 2: Is/It<1, Stable, Exponential Model 
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Figure (21) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 2: Is/It<1, Stable, Exponential Model 
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Figure (23) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 2: Is/It<1, Stable, Exponential Model 
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Figure (25) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 2:Is/It<1, Stable, Exponential Model 
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CASE 3: 4s <1 , Unstable, Exponential Model 


t 


System Parameters Initial Conditions 


Platform Rotor Platform Rotor 
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Table (4) Case 3 Parameters 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Figure (27) Total Energy and Percent Difference Angular Momentum Versus Time 


Case 3:Is/It<1, Unstable, Exponential Model 
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Figure (28) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Figure (29) Total Energy and Rotor Core Energy Versus Time 
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Figure (30) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Figure (31) Total Energy and Platform Core Energy Versus Time 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Figure (34) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 3: Is/It<1, Unstable, Exponential Model 
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Figure (36) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 3: Is/It<1, Unstable, Exponential Model 
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CASE 4: : >1 , Unstable, Exponential Model 
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System Parameters Initial Conditions 


Platform Rotor Platform Rotor 
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c= 1.0 sec C= 1.0 sec 


Ee final = 1247.1J Eefinat = 1552.5] 
Ec bic 


fe 03m = stot = 1.025 
Te total 





= positive 


= ~,0005 s-! 





Table (5) Case 4 Parameters 
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Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (38) Total Energy and Percent Difference Angular Momentum Versus Time 


Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (39) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Is/It>1, Unstable, Exponential Model 
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| Figure (40) Total Energy and Rotor Core Energy Versus Time 


Case 4: Is/It>1, Unstable, Exponential Mode’ 


core energy (J) 





Figure (41) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (43) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (44) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 


Case 4:Is/It>1, Unstable, Exponential Model 
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Figure (45) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (47) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 4: Is/It>1, Unstable, Exponential Model 
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Figure (48) Core Energy Versus Nutation Angle 
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CASES: <1 , Stable, Verhulst Model 


t 


System Parameters Initial Conditions 
Platform Rotor Platform 


0.0 
hh’ = 1600 kgm? = 1000kgm? - 


I3’ = 1500 kg m2 1200 kg m2 
M’ = 1000.0 kg 700.0 kg 


m= 20.0 kg 1.0 kg 
1.0 m 1.0 m 


m 
0.0 © 


oe 5 rad 
@3 = 7.27 x 10° cer 





rad 
0.10 So 


Core Energy Parameters 


Platform Rotor 


N N 
= k= 1.0 = 


kg kg 
WM ee 1.0 Se 
rie — See I) E- initict = “Sas 


Ecfinal’ = 0.1223  Ecpinat = 3170.93 
Ec’ 

4 
r = -.0002 s"! -.0002 s-! 


Estoial = 01.907 
Ts total 





. E ; 
= negative negative 





Table (6) Case 5 Parameters 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (49) Total Energy and Percent Difference Angular Momentum Versus Time 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (50) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (51) Total Energy and Rotor Core Energy Versus Time 





Case 5:Is/It<1, Stable, Verhulst Model 
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Figure (52) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (53) Total Energy and Platform Core Energy Versus Time 





Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (54) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (56) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (58) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 5: Is/It<1, Stable, Verhulst Model 
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Figure (59) Core Energy Versus Nutation Angle 
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CASE 6: i <1 , Unstable, Verhulst Model 


t 





System Parameters Initial Conditions 
Platform Rotor Platform Rotor 

0.0 m 

1’ = 1600 kgm? I,= 1000kgm2 a 
0.0 

3’ = 1500 kgm? Iz= 1200kgm? s 

, 5 rad 

f’= 1000.0 kg M= 700.0 kg 7.27 x 10™ Se 

n= 1.0 kg m= 1.0 kg ra 
0.10 tad 

7= 10m a= 1.0 m 

Y= 1.0 ot k= 1.0 a Core Energy Parameters 

. kg - kg Platform Rotor 

c= 1.0 sec c= 1.0 sec 


es a Pain 3061.6 J 


L= 10m ‘stol= 0.897 1161.03 — Ecfinat = 1489.9 


I t total 





. E x 
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-1x 10% s-! ~1x 10% 5-1 





Table (7) Case 6 Parameters 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (60) Total Energy and Percent Difference Angular Momentum Versus Time 


Case 6: Is/It<], Unstable, Verhulst Model 
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Figure (61) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (62) Total Energy and Rotor Core Energy Versus Time 


| Case 6:Is/I:- Unstable, Verhulst Model 
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Figure (63) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (65) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (67) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (68) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 


Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (69) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 6: Is/It<1, Unstable, Verhulst Model 
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Figure (70) Core Energy Versus Nutation Angle 
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CASE 7: 7 >1 , Unstable, Verhulst Model 
f 
















System Parameters Initial Conditions 
Platform Rotor Platform Rotor 
0.0 
I1,= 1000 kgm? m 
0.90 mM 
z= 1200 kgm? S 
5 rad 
M= 700.0 kg 7.27 x 10° a 
m= 1.0 kg 
0.10 fad 
a= 1.0 m 
k= 1.0 = Core Energy Parameters 
_ kg Platform Rotor 
C= 1.0 sec ee vee a On nce ee 
Ec = 13.2 J E this 3059.7. 
Sta = 1,025 Ec final’ = 1247.13 Ecfinat = (1552.5. 
t total — : 
a positve te positive 
NV’ A 


~1.1x 10% s-! r = -1.1x10~ 









Table (8) Case 7 Parameters 
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Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (71) ‘Total Energy and Percent Difference Angular Momentum Versus Time 


Case 7:Is/It>1, Unstable, Verhulst Model 
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Figure (72) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (73) Total Energy and Rotor Core Energy Versus Time 


Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (74) Total Energy and Rotor Core Energy - First 100 Seconds 
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Case 7: Is/It>1, Unstable, Verhulst Model 


----- platform core energy 


ew ne 


0 1000 4000 5000 


we 
ot od 
poww eon nnen rr” 


Figure (75) Total Energy and Platform Core Energy Versus Time 


Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (76) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (77) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 


Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (78) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (79) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 


Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (80) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 7: Is/It>1, Unstable, Verhulst Model 
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Figure (81) Core Energy Versus Nutation Angie 
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APPENDIX B COMPUTER PROGRAM CODE 


a) 
4 


anna t ewe eon ee eee se declarations RaGARARREHRA RO RE AH AR EADS / 


#@include <math.h> 
#include <malloc.h> 
finclude <stdio.h> 


édefine MAXSTP i500c0 
tdefine TINY 1.0e-S¢ 


#define PGROW -C.20 
édefine PSHRNK -0.23 


/* odeint °*/ 
/* adeint *f 


is brqew*s 
ad oagie i! 


#define FCOR 0.06665666666666666666666666666667 /* 1.0/15.0 krqc*/ 


#define SAFETY 0.9 
#define ERRCON 6.0e-4 


*xx=0; Ties 
kount=0; libs 
**ype0, dxsave0; oa 


float **y=-0, 
int kmax=0, 
float *xp=0, 


[/etteretenateteorere Grror function 
void nrerror(error_ ext) 
char error _texti]; 


{ 


void exit(); 


7° kKeqe-*/ 
Leaakrac << / 


defining declaration rkdumb ‘*/ 


defining declaration odeint */ 
defining declaration odeint */ 


Reehekthe he theaeeka han t / 


fprintf(stderr, "Numerical Recipes run-time error...\n"); 


fprintf(stderr, "ss\n",error_ text): 
fprintf(stderr," .. 
exit(l): 


.now exiting to system. 


apn iba) es 


[/eat*taeeeehaaaeeoonteran vector function ReARRAT HARRAH KRHA HERRERO / 


float *vectorr(nl,nhj 
mene ni, nh: 


{ 
float *v:; 


ve (float *)malloc({ (unsigned) 


(nh-nl+1)*sizeof (float)): 


4£ (!'v) nrerror("aliccation failure in vector()"); 


return v-nil: 


) 


« * = 
ii Rea akRneE RRA KRAREeTRE WS matrix function ARkeeEHHAHAEERHERA REE E AEA SE OO / 


float 
int nrl, 
{ 

mnt i: 
float 


mem, nel, -.c2A: 


Bem: 


me (float **) malles{ (unsigned) 


(nrh-nrl+t+l)*sizeof (float*)); 


1£ (!m) nrerror("aliccation failure 1 in matrix()"); 


m-= nril; 
for (ienrl;i<enrh:i+-) 


{ 


m{ij = (float *) malloc( (unsigned) 
mrerror ("allocation failure 2 in matrix()"): 


wee im[i;} 
mii} 
) 
return m; 
j 


[/teeawanennenenetenee 
void free vector (v,ni,ah) 
Zloat *v:; 

mac nl, nh: 


(nch-ncl+l)*sizeof (float)): 


£ree vector function SHARARAHREARKREREHREEASE HOS / 
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free((char*) (vtnlige 


) 


[asset r ee eaee ne ne eee t ees free matrix function ee Ae ee ee ef 
void free_matrix(m,nrl,nrh,ncl,nch) 

float **m; 

int nrl, narh, neieench; 


{ 
Pntests 


for (ienrh i >=nxl: ia free ((char*) (m[i)4#ncl)); 
free((char*) (m+tnrl)); 


} 


[/ReKkawetetkeeeAsanenes rk4 function HAHARTHERHERHRABRAR HE EASE / 


void rk4{y,dydx,n,x,h, yout, derivs) 


float yi). dydx({). x, Ho youtl)-: 
void (*derivs) ():; 
inten: 


{ 

inti: 

float xh, hh, hb; *dym, “dyn cave yc. VECtCOrri{): 
void fee vector GE 


dym = vectorr(l1,n);: 
dyn = vectorr(1,n): 
dyt = vectorr(1,n): 
yt = vectorr(l,n); 
hic= hw 07 5. 
h6 = h/6.0; 
xh = xthh; 
for (iplsisen;it++) yt(i) = yli)t+thh*dydx{iJ; /* first step a, 
(*derivs) (xh, yt, dyt,dydx); 
for (iel;i<en;i++) ytli) = yli)J+hh*dyt[iJ-; 
(*derivs) (xh, yt, cym, dyt); 
for (181; i¢c=n;1+7) 
{ 
yt(i) = ylajehtdymiae 
dyn[(i) = ecyt(i)+dym(i); 
) 
(*derivs) (x+h, yt, c;t, dym) - 
for (ielsi<en-;it+) yout[i]) = y[i)+h6* (dydx[iJ+dyt[i]+2.0*dyn[ij): 
free_vector(yt,1,n); 
free.vecctor (dyt, 1 n}: 
free vector(dym, 17 n); 
free veccor(dyn, 1, n)- 


) 


fetteaeewnetetetanetreeese rkdumb function HHEETHHHRARHREHHHEREHA EAE HO od / 


void rkdumb(vstart,avar, xl, x2,nstep,derivs) 
int nvar,nstep; 

float vstart |), x17x2: 

void (*derivs) (); 

{ 

int -i5 7k: 

float x, fh: 

float *v,.-*vout,> *cv *vectorm a. 

void rk4(), mrerrox(), free_vector(): 


v = vectorr(l,nvar): 


vout = vectorr(l,nvar}): 
dv = vectorr(l,nvar): 
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for (i#l;i<envar; i++) 


v[i) = vstart[i): 
vues s) = vlil: 
) 
xx{1}) = xl; 
x = xl; 
h = (x2-x1l)/nstep; 
for (kel; k<enstep; k++) 
; { 
(*derivs) (x,v,dv); 
rk4(v,dv,nvar,x,h,vout,derivs):; 
if (xth *= x) nrerror("Step size too small in routine RKDUIRN): 
x += h; 
xx(k+l] = x: 
for (lel; i<envar; i++) 
( 
v[i) = vout[i); 
wirik+i) = v (ij; 
) 
) 
free_vector (dv,1,nvar); 
free vector (vout,1,nvar); 
free vector (v,1,nvar); 


) 


ier we eee ae an aane st rkqc function RAREKARKERAKAAARKARR RRA RHA ARERR / 


void rkqc(y, dydx,n,x,htry,eps, yscal, hdid, hnext, derivs) 


muoaney|), Gydx[), *x, htry, eps, yscal[]), *hdid, *hnext; 
void (*derivs) (); 
mat nN; 


( 

nt. i; 

float xsav, hh, h, temp, errmax: 

float *dysav, “*ysav, *ytemp, *vectorr(); 
void rk4(), nrerror(), free _vector(); 


dysav = vectorr(l1,n):; 
ysav = vectorr(1,n): 
ytemp = vectorr(i,n): 
msavy = (*x): 


for (iel:i<en; i++) 
{ 
meav(i) = y[3)}: 
Cysav[i) = dydx[i); 
) 


he htry; 
more ; : ) 
( 
Pes 0.5*h: 
rx4(ysav,dysav,n,xsav,hh, ytemp, derivs);: 
*x = xsav+hh; 
(*derivs) (*x, ytemp, dydx); 
rk4 (ytemp, dydx,n,*x,hh,y,derivs); 
*x = xsavth; 
if (*x = xsav) nrerror("Step size too small in routine RKOC"); 
rk4(ysav,dysav,n,xsav,h,ytemp,derivs) ; 
errmax = 0.0; 
for (isl; i<en; i++) 
{ 
ytemp[i]) = y[iJ-ytemp[i]J-: 
temp = fabs (ytemp[iJ]/yscal[i)); 
if (exrmax < temp) errmax = temp; 
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) 

errmax /= eps; 

if (errmax <= 1.0) 
{ 
*hdid = h; 


*hnext = (errmax > ERRCON ? SAFETY*h*’exp(FGPOW! loa (errmis)) 


break; 
) 
h = SAFETY*h*texp (PSHRNK*1log (errmax) ); 
) 
for (iel;i<en:it+) y{i}) += ytemp[i]) *FCOR: 
free_vector(ytemp,1,n); 
free vector (dysav,i,n); 
free vector (ysav, lin); 


) 


[tee nenaneeneneneeoes odeint function eeeeeaeeeeeeeh ae we eases / 


void odeint (ystart,avar,x1l,x2,eps,hi,hmin, nok, nbad, derivs, rkqec) 


float ystart{], xl, x2, eps, hl, hmin; 
int nvar, *nok, *nbad; 
void(*derivs) (); 

Voud (*rkGe) i): 


{ 

Int ast pea. 

float xsav, x, hnext, hdid, h:; 

float *yscal, *y, *dydx, *veccoun«): 
void nrerror(), free vector (); 


yscal = vectorr(l1,nvar); 

y @ vectorr(l,nvar): 

dydx = vectorr(l,nvar); 

x = xl; 

h = (x2 > xl) ? fabs(hl) : ~fabs(hl); 
*nok = (*nbad) = kount = 0; 

for (iel;i<envar;:i++) yf{i}) = ystart{(i]; 
if (kmax > 0) xsav = x~dxsav*2.0; 

for (nstpel;nstp<=MAXSTP; nstp++) 


(*derivs) (%,y,dydx); 


for (iel:i<envar:i++) yscal(i) = fabs(y(i))+fabs (dydx[i)*h) ¢7Tiur: 


if (kmax > C) 
{ 
if ‘fabs(x-xsav) > fabs (dxsav) ) 
{ 
if (kount < kmax-i) 
{ 


xp({++kount) = x; 


for (ielsis<=nvar:i++) ypl{i)ikount)—-; ia 


xSav = x; 
) 
) 
) 
if ((x+h-x2)° (x+h~-xl) > 0.0) h = x2-»; 
(*rkqc) (y,d;¢cx,nvar,&x,h,eps, yscal, &hdid, &hnext,dcrivs): 
if (hdid ©== &) ++(*nok); else ++(*nbad); 
if ((x-x2)* (x2~-x1) >= 0.0) 
{ 
for ‘iel:i<envar; i++) 
ystart(ij] = y{ij: 
if (xmax) 


{ 


xp(++kount)©x; 


for (iel;i<envargit+) Pyota) kounejecey (1). 
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4 


fit tp ye 


) 
free vector (dydx,1l,nvar); 
free vector (y,1,nvar); 
free vector (yscal,1l,nvar); 
return: 


} 


if (fabs(hnext) <= hmin) nrerror("Step size too small in QnDEINE™): 
h = hnext: 
: 


nrerror("Too many steps in routine ODEINT"); 


) 
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[/AeARHRAHKEARARAATAS ES declarations 


finclude 


AaRAERARAEHRHARAHRRHE AR HHA HARA HE HEEB DEY 


a dd.) lh 


#fdefine N ? 
#define MAXARRAY 100000 


#fdefine SQ(x) 


((x)* (x)) 


int kount, kmax;: 

figat **vp, *xp; 

float dxsav, dxmin, dxlimit: 

float I1=-0.0, I2°0.0, I3=0.0, M=0.0, m=0.0, mb=0.0, a=0.0, k-90.0, c-0.%: 
float Ilp=0.0,I2p=-0.0,I13p=0.0,Mp=-0.0,mp=0.0,mbp=0.0,ap-0.0,kp-O.0,ep 0.0; 
float L=0.0, sigma=0.0, MT=-0.0, v=0.0, rho=0.0, rhop-0.0, L1-0.0, 1.2-0.0: 
float A-0.0, C=#0.0, J3p=-0.0, w3p=0.0; 


[MERAH EERE EERE EOE EQUATIONS OF 


MOTION WAAR HARHEHA ERE HEE EEE EAE D, 


void derivs(x,y,dydx) 


float x, y{J, dydx[]:; 

{ 

float zeta=-0.0, dydxzeta=0.0; 

float Al=0.0, A2=0.0, A3=0.0, A4=0.0, AS5=0.0, A6=0.0, AT-C.0, AR-O.0: 
float A9=0.0, Al0=#0.0, Al1=#0.0, Al2=0.0, Al13=0.0, Al4=0.0, A15-0.0, AI6-9.0 
float Al7=-0.0, Al8=0.0, A19=0.0, A20=0.0, A21=0.0, A222#0.0, A23B-N.0, ATA-O. 
float A25=0.0, A26=0.0, A27=#0.0, A28=0.0; 

float Bl=0.0, B2=0.0, B3=0.0, 8B4=-0.0, BS5=0.0, BG6="0.0, B7-0.9, FAR-0. 
float B9=0.0, Bl0=0.0, Bll=0.0, B1l2=0.0, B13=0.0, B14=0.0, £215-0.0, fR16-0.0 
float Bl7=-0.0, B1l8=0C.0, B19=0.0, B20=0.0, B21=0.0, B22-0.0, B23-0.9, BGO. 
float B25-0.0, B26=0.0; 

float C1=0.0, C2=6.0, C3=0.0, C4=0.0, C5#0.0, Cc6=-0.0, C7-N.0, C#-4,0; 
float C9=0.0, C10-0.0, C1l=0.0; 

float Dle-0.0, D2-0.0, D3=0.0, 0D4"0.0, D5“0.0, D6=0.0, D7-0.0, NS-S.8 
float £1=0.0, E2-C.0, E3°0.0, E4"0.0, £5=0.0, E6=0.0, E7-0.9, E&-0 
float E9=0.0, £10=-0.0; 

float Z1=0.0, Z2-0.0, 2320.0, 24="0.0, Fl=0.0, F2=0.0; 

dydx{4J = y{5]; 

dydx{6}] = y[7)-; 


zeta © rho*y(4)+rhop*y(6); 


dydxzeta = rhoty([S]¢rhop‘yi 7]; 


Al = -(A-6) viz eal); A2 = J3p*sigmaty(2); 

A3 = 2°*MT*zeta*dydxzeta*y([1l]}; A4 = MT*SQ (zeta); 

AS = -MITSOtzeta) syi2 ity isle AG = -2*m*(zeta+L2)*y{[4]: 

AT © 2*m* (zeta+L2) *v {4} *y[2]} *y [3]: AB = -2*m* (zetatL2)*y[5)*;[1]: 

WS = =2° my 4) covdxcetasy pL) AlO = 2*mty{4)*y(S5)*y[1]): 

All = m*SQ(y[4)J); A125 = =m" SOC yal) ya ele 

Al3 © -m*ty{[4)*a; AG = -m* yl 4itaty (ey i 

AlS = -2*mp*(zeta—=b1)*y (6); Al6 = 2*mp* (zeta-Ll)*y([6)‘y 12) ‘y[3): 
Al? = =<2*°mp*(zeta-c1)-y (7 lovin A188 = -2*mp*y(6)*dydxzeta’y [1]: 

Al9 = 2*mp*y16)*vi7)}>y11)1- A2Q = mp*SQ(y[6)): 

A2l = -mp*SQ(y[6)) *y[2])*y{3)-; A2Z2 = =-mp*y([6] *ap*cos (siqza’s): 

A23 = -mp*y[6)*ap*cos (sigma*x) *y(1])*y[2]; AZ4 = mp*’ap Si -i(siGqira “2 ie 
A25 = mp*ap*sin(sigma*x) *y[{6)* (SQ(y [3] +sigma)-SQ(y[2])); 

A26 = AtA4+A6+A11+4A15+A20; A27 = A13+A22; 

A228 = A1+A2+A3+AS+AT+AB + AS4A1LO+ALZ+ AL 4+ Al G+Al T4HALB+ALVIAZIAAZ3AA2S: 

Bl. © = (CHA) yt ew Si: B2 © -J3p*sigmaty{ i); 

B3 = 2*MT*zeta*dydxzeta*y(2]j; B4 = MT*SQ(zeta): 
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BS = MT*SQ(zeta) *y(1j*y(3):; B6 = -2*m* (zetatL2) *;([4;: 
B7 = ~-2*m* (zetat+L2) *y[(4) *y(1)*y(3); BB = -2*m* (zeta+L2) *y(Sicyl7l: 
B9 = -2*m*y[4) *dydxzeta*y [2]; B10 = 2*mty[4)*y([5)*yi?7]- 
Bll = m*SQ(y[4)): Bi2 = m*SQOty[4)) * yt) *yi*:: 
B13 = -m*a; Bl4 = -mtaty (4) *(SOl(y(3}) -SOly ll )))- 
B15 = -2*mp* (zeta-Ll) *y(6]); B16 = -2*mp*(zeta-L])‘y[& ‘yllltyt?l: 
B17 = -2*mp* (zeta-Ll) *y (7) *y(2); B18 = 2*mp*y(6)*yl2)* (yl 7. -dydssetad: 
B19 = mp*SQ(y[6]): B20 = mp*SQ(y[6J)‘yli)*yl 7): 
B21 = -mp*ap*cos(sigma*x) ; B23 = -mp*ap*sin(siagma’~) ‘y [6]: 
B22 = -mp*tap*cos(sigma*x) * (SQ(y[3]+sigma) -SQ(y[1])) *y{6]; 
B24 = mp*ap*sin (sigma*x) *y(6] *y[1])*y(2); B25 = AtH4 AG VIII wie: 
B26 = B1+B2+B3+B5+B7+B8+B9+B10+B12+B14+B16+B1l7+B18+B20+B22+h24: 
Cl = -2*m*taty(5) *y(1)- C2 = -mtaty([4)-; 
C3 = mtaty(4)*y(2) *y(3); C4 = -2*mptap*sin(sigma’=jty[7}‘ylel. 
C5- = -mp*ap*sin(sigma*x) *y[6); C6 = -mp*tap*sin(sigma*x)*y[Cl}fyllliyl '-: 
C7 = -2*mp*ap*cos(sigma*x)*y(7)*y[1}; C8 = -mp*tap*cos(sigma*s) *y [6]: 
C9 = mp*ap*cos (sigma*x) *y(6)*y[2)*y(3]; ; 
C10 = C2+C8; Cll = C1+C34+C4+C6+C7+C9; 
D1 = m*(l-rho); D2 = -mp*trho; 
D3 = ~-m*a; D4 = mtaty([i)*y(3]-: 
DS = -m* (SQ(y[1))+SQ(y(2]))* (y(4] * (l-rho)-L2-rhop*y[6)); ews ylal 
D7 = k*y[4); D8 = D4+DS+D6+D7; 
El = -m*rhop; E2 = mp* (l-rhop): 
E3 = mp*ap*sin(sigma*x); E4 = mp*ap*sin(sigma*x) *y[2]* (yi 3) 62° sien): 
ES = -mp*ap*cos (sigma*»x) ; E6 = mptap*cos (sigma*x) *y[1l)*(y{[3)'2¢*siari): 
E7 = -mp* (SQ(y(1))+SQ(y[2]))*(y[6]* (l-rhop) +Li-rho*y([4]); FR = ep yl id: 
E9 = kp*y([6]; E10 = E4+E6+E7+ES+E9: 
enetA2 7 *B1 3) / (AZ6*B23) -E1/E3)/ ((A27*B25) / (A26*B23) +E5/E3) ; 
Z2 = (-(A27*B21) /(A26*B23) +A24/A26-E2/E3) / ((A27*B25) / (A26*B23) 45/53) ; 
Fl = (-(A27*B26) / (A26*B23) +A28/A26-E10/E3)/( (A27*B25) / (A26°B23) 425/53): 
ZS @ ((C*B13) / (C10*B23) - (A27*B13) / (A26*B23))/ 

( (A27*B25) / (A26*B23) - (C*B25) /(C10°R23)+C5/C10): 
Z4 = ((C*B21)/ (C10*B23) -(A27*B21) / (A26*B23) +A24/A26) / 


(127 *B25) 7 (AZ6*B23)-(C*B25)7 (C10°R2 3) 4+CS/CIC): 
F2 = ((C*B26)/ (C10*323)-C11/C10=(A27*B26) / (A26*B23) +A28/A26) / 
(ONZ TR B2 5) (N26 B23) -(C*B25) 7 (Ci0*v? 3) 1Ch/C10): 


dydx(5) = ((-F2-D8/D3)/ (24+D2/D3) + (F1+D8/D3) / (22+D2/D3) )/ 
((-Z1=-D1/D3) / (Z2+D2/D3)4+(Z34D1/N3) / (7 41D2/D Ay): 


dydx(7) = ((-F2-D8/D3)/ (23+D1/D3) + (F1+D8/D3) / (21+D1/D3) )/ 
((-22-D2/D3) / (Z1+D1/D3) + (244D2/D3) / (7341 /I2)): 


dydx(3) = (dydx{5)* (C5*B13)/(C*B25)+ 
dydx [7] *((C10*A24) / (C*A26) + (C5*B21) / (C*B25) ) + 
(C10*A2B) / (C*A26) + (C5*B26) / (C*B25) -C11/C) / 
(1- (C10*A27) / (C*A26) - (C5*B23) / (C*B25)): 


dydx(1} = dydx[3) * (-A27/A26) +dydx (7) *(-A24/A26) + (-A28/A26) ; 

dydx[(2) = dydx{3) * (-B23/B25) +dydx(5] * (-B13/B25) +dydx (7) * (-B21/B25)-P26/R25: 
) 

[tte tenatateeetaeatavenne MAIN PROGRAM Ree ee See ee A SA ER ARO FP AER ets 8 8 ff 
main () 


{ 

ant eegeenoad, nok: 

enc 11=1000, 12=20C9, 13=3000, 14=4000, 15=5000, 16=6000, 17-7098: 

ant 18=8000, 19=9000, 1100-10000; 

s20at eps, hi, hmin, dxsavd, xl, x2, xscale, *ystart; 

float H1"0.0, H2=0.0, H3=0.0, Hlp=0.0, H2p=0.0, H3p<0.0, zem-0.0, A: gumen: : 


ize 


float T1-0.0, T2-0.0, T3=0.0, T4=0.0, T5=0.0, 7620.0, 1) {0B 24. 
float T9=-0.0, T10°0.0, T11=#0.0, T12=20.0, T13=0.0: 

float zemrigid, It, Itp, Ittotal,@is, Ispyeistotal, wn, fanaa ia, 
float *ke, *kedel, *etotal, *etotaldel, ‘*h, *hdel, *theta; 

float *Ec, *Ecdel, *Ecp, *Ecpdel, *Echype, *Echypedel, ‘“*Ecphypo, ‘*Fephypede 1; 
fioat Q, Op, *eta, *etap, Eemin, ecfactor, ecpfactor, hsq, "betina! et =I sie 
float signec, signecp, signg, signqp, arg, argp, Ecf, Ecfp: 
float *Ecexp, *Ecpexp, *Ecexpdel, *Ecpexpdel, ‘*etah, *etahp: 


ystart = vectorr(1,N); Ec = vectorr (i, ress APPAT): 
ke = vectorr(1,MAXARRAY) ; Ecpdel = vectorr(],MAZARRAY) ; 
kedel = vectorr (1,MAXARRAY) ; Echype = vectorr (1, 1AZARPAY); 
etotal = vectorr (1,MAXARRAY) ; Echypedel = vectorr(}, (ANAFPRAT) ; 
etotaldel = vectorr(1,MAXARRAY) ; Ecphype = vectorr ()] i sARPAYT): 
h = vectorr (1,MAXARRAY) ; Ecphypedel = vectorr(1],tin“Aankar) ; 
theta = vectorr(1,MAXARRAY) ; eta = vectorr( (1, nsARRAY): 
hdel = vectorr(1,MAXARRAY) ; etap = vectorr(1],iA2zARRAY); 
Ecexp = vectorr(1,MAXARRAY) : Ecexpdel - vectorr (1, AZ ARBAY) ; 
Ecpexp # vectorr(1,MAXARRAY) ; Ecpexpdel = vectorr(],°:neAEPAY); 

etah = vectorr(1,:722ARRAT);: 

etahp = vectorr(liiecARRAT): 


scant ("SEVESEVLLELLREVLVEVERLVELERLRELLRESLCVEVLRCVELSELESCSLSCAL 2070", 
&x2,611,613,6M, im, 6a, &k,6C,6L,&I1p, &13p, &Mp, &mp, ap, &k pp, itp, 
éwop,éystart(l),éystart(2),éystart(3),éystart[4),&ystar[5)., 
éystart[(6),éystart[(7),éeps, &6ecfactor, Secpfactor,ésEcf,fe-ir): 


x1=-0.0001; h1i-0.00001; hmine1.0¢e-11; kmax=-100000; 
mb = m; mbp = mp: I2-= 11; I2pe-) lip: 
sigma = w3p-ystart([3); 


/**** critical damping - input damping = 1000.0, makes it crit:771 ***?*/ 
if (c¢ ©= 1000.0) ¢ = 1.5*sqrt (4*m*k); 
if (cp == 1000.0) cp = 1.5*sqrt (4*mp*kp); 


/*** establish time interval for saving data ***/ 
dxmin - 1.0e-4; 
Gxlimit = (x2-x1)/1750.0; 
dxsav = (x2<20.0) ? dxmin : dxlimit; 


Ma, = M+Mp+4.0*mb+4.0*mbp; 

v = (Mp+4.0*mbp) /MT; 

A © Il+Ilp+2.0*m5*SQ (a) +2.0*mbp*SQ (ap) + (Mpt+4.0*mbp) * (1.0-v) °77 (i): 
c ~ I3+I3p+4.0*m5*SQ (a) +4.0*mbp*SQ (ap): 

J3p © I3p+4.0*mbp*SQ (ap); 

rho = m/MT; 


rhop = mp/MT; 
Ll = L* (M+4.0*mb) /MT: 
L2 = L*(Mp+470*mbo) /MT; 


zemrigid = ((Mp+4.0*mbp) *L) /MT: 


ie = I1+M*SQ(zemrigid) +mb* (4.0*SQ(zemrigid) +2.0*SQ(a)): 

itp ~ lip+Mp*SQ(L-zemrigid) +mbp* (4.0°SQ(L-zemrigid)+2.0°S¢ ‘=:)}: 
attotal = [tv ith- 

Is = 13+4.0*md*SQ (a); 

Isp = I3p+4. 0*md>p*SOtap)- 

Istotal e- Is+tIsp; 


odeint (ystart,N,x1l,x2,eps,hl, hmin, 6nok, 6nbad, derivs,rkqe): 
for (jel; J<@kount; j++) 

( 

zem = ((Mp+4.0*mbp) *L+m*yp([4} [j]+mp*yp[6) [j)) /MT:; 


HL = (11+M*SQ(zem) +m* (2.0*SQ (a) +4.0*SQ(zem) +SO(ypl4) (55: 
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H2 


B3 


Klip 


H2p 


H3p 


-2.0*zem*yp([4) (3))) *yp(1) [j)-m*a*typl4) (35) typl3) (3): 


(I24+M*SQ(zem) +m* (2.0*SQ(a) +4.0*SQ(zem) +SQ(yp[4) (5)) . 
-2.0*zem*yp (4) (3))) *ypl2) (j)]-m*ta*ypi5) (5): 


-mta*yp (4) (4) *yp{1) (4) +(13+4.0*m*SQ(a)) *yp[3) [5]; 


(Ilp+Mp*SQ(L-zem) +mp* (2.0*SQ(ap) +4.0*SQ(L-zem) +SQ(ypi6) [5]) 4 
2.0* (L-zem) *yp[6) (3))) *yp(1) (4) -mp*ap*cos (sigma*xp|[))) 
yp(6) (3) *yp[3) (4) -mp*ap*sigma*cos (sigma*xp[))) * 
yp[6) [j)+mp*ap*sin(sigma*xp(4)) *yp{7) (Jj): 


(I2p+Mp*SQ(L-zem) +mp* (2.0*SQ(ap) +4.0*SQ(L-zem) +SQ(yp[6] [3))! 
2.0* (L~zem) *yp[6) [j))) *yp(2) (3) -mp*ap*sin(sigma*»p| 5] )* 
yp(6) (3) *yp[3) [j) -mp*ap*sigma*tsin (sigma*xp|[j)) * 
yp(6) (j)-mp*ap*cos (sigma*xp[j)) *yp[7) (3): 


-mp*ap*cos (sigma*xp[j)) *yp(6) (3) *ypl1) (3)}-mp*ap* 
sin(sigma*xp([4j))*yp(6) (4) *yp(2) (3) +(13p+4 .0*mp*S0O (ap) ) * 
(yp{3) (4) +sigma) ; 


h[(3) = sqrt (SQ(H1+H1p) +SQ(H2+H2p) +SQ(H3+H3p) ); 


argument = (H3+H3p)/h([j): 


if (argument>0.99999999) argument = 1.0; 


theta[j) = 57.2957795131*acos (argument) : 


hdel(3) = h[j)-hl1): 


a1 


T2 


T3 


T4 


iD 


T6 


ne 


T8 
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T10 


pal 


m2 


0.5*(I1*SQ(yp (1) [j))+I2*SQ(ypl2) (3))+13*SQlypl3)[5))): 


Oro" (Tip*SO(vold) (7) 42 2p*sO(ypl2Z).15))+I3p* 
SOtyp(3))\( 314 sikema)) 


0.5*m* (SQ(yp[5) [3))+(SQlypl1) (3)) +SQ(yp[2)(5))) *SQ(yp[4) Li)) 
=—An0caevP (2) (9) *ypis)(j)*yel4) ii}. 


0.5*mp* (SO(yp[7) (3) )+SQl(ypl6) (5)) * (SQ(yp[1)(53))+SQ(ypl2} 05))) 
+2.0*ap*sin (sigma*xp[j)) *yp(2) (3) *yp(3] (i) *ypl6) (3) 
-2.0*ap*cos (sigma*xp[j))*yp(1) (3) *yp(3) (3) *yplS) (531): 

emeatyp[2) (3) *yp iS): 

“mp* (-~ap*yp[7) [j) *yp(2) [j) *cos(sigma*xp[)j)) 
+ap*yp(7) (53) *ypl1) (3) *sin(sigma*xp[j)) 
~ap*sigma*typ([1) (3) *yp[6) [j}] *cos (sigma*xp[j)) 
~ap*sigma*yp([2) [j)*yp[6) [j) *sin(sigma*xp[j])):; 

(0. 5*SQ(m*yp([5) (j)+mp*yp([7) (3)))/MT: 


(Omo* (SOlyp(1)[j))+SOlyp{2) 13))) 
*SQ(m*yp[4]) [j) 4mp*typ(6) (5) )) /uT: 


0.5* (SQ(yp[1) (4))+SQ(yp[2) [5))) * (M+4.0*m) *SQ(L) * 
SQ((Mp+4.0*mp) /MT): 


0.5*(SQ(yp(1) (3))+SQlyp(2)13)))* (Mp+4.0*mp) *SQ(L) * 
SQ((M+4.0*m) /HT); 


(m+mp) * (-yp(5) (3))* (m*yp[5) [j)+mp*yp(7) (3)) /MT: 


-m*yp (4) (4) * (SQ(yp[1) [j))+SQ(yp(2) (5)))* 
(m*yp[4) (3)+mp*typ[6) (3) + (Mpt+4.0*mp) *1,) /tiT: 
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T13 = -mptyp(6) (3) * (SQ(ypl1)(5)) *SQ(ypl2) (31))* | 
(meypo(4) (Jj) 4mptyp(6) (jy- (tie4 Oem) Ob) ZU: 


T14T24T9tT4t+THt+TOFT 7 +TS8ITOIATIOFTI1IIATIA (TI i: 


ke[j] 


kedel[j]) ke(j)-ke[1); 


etotal[}) = ke{j)+0.5*k*SQ(yp(4](j)) +0.S*kp*sSQlypl6} (i)); 
etotaldel(j) = etotal[j]-etotal([i]- 


Ec[ 3) = 0.5*(Ittotal*(SO(ypi1) [jj])4+SO(yptz) lee . 
+Istotal*Sso(ypf i} ild) 


Ecdel[j) = Ec{j)-Ec(1]> 


Ecp[j) <= 0.5* (1ttotal* (SQlypi1)}19))4S0lyel7 tae | 
+lstotal*sSsQ(yp(3) (i) 's1ama)): 


Ecpdel[{j) = Ecp(j)-Ecpl1]: 
} 


wn = (Is*yp{3)(1)+Isp* (yp(3) (1)+sigma)) /Ittotal; 

lambda = wn-yp([3) (1); 

lambdap = wn-(yp([3) [1])+sigma) ; 

hsq = SQ(Ittotal) *(SQ(yp(1) [1])+SQlyp(2)(1))) 
+SQ(Istyp(3) (lJ) +Isp*(ypl3) tl) 'siqma)): 


if (theta[kount)-theta[l] < 0.0) /**** stable condition ****/ 
{ 
if (lambda < 0.0) { signec = 1.0; signg = 1.0: ] 
else ( signec = -1.0; signgq = -1.0; } 
1£ (lambdap < 0.0) ( signecp = 1.0; signgp - 1.0; } 
else { Signecp. = -1.0°-.signep = -)..0 
] 

else /**** UNnStable Conditions. 7 
{ 
if (lambda < 0.0) ( signec = -1.0; signg = 1.0; 3 
else ( signec = 1.0; signg = -1.0; ] 
if (lambdap < 0.0) { signecp = -1.0; signqp = 1.0: ] 
else ( signecp = 1.0; signqp = -1.0: ) 


) 


Ecfinal = (Ecf w= 0.0) ? Ec(kount) : Ecf; 
Ecpfinal = (Ecfp "= 0.0) ? Ecp[(kount]) : Ecfp; 


for (jel: j<ekount; 3++) 
( 
Echypel[j) © Ec(1) *Ecfinal/ (Ec(i)+(Ecfinal-Ec{1)) * 
exp(ecfactor*Ecfinal’xp[j]})): 
Echypedel(j) = Echype[j) -Echype(1): 


Ecphype [34] - Ecp([(1l)*Ecpfinal/ (Ecp(1)+(Ecpfinal-Eep([)))>- 
exp(ecpfactor’*Ecp:finaltezplj))): 


Ecphypedel[j) = Ecphype[j]-Ecphype(1); 

Q = 2.0*Echype (jj) *Istotal*(1.0-(Istotal/Ittota}) ) 
-~hsq* (Istotal/Ittotal)*(1.0-(Istotal/Ittotal)) 
+(Istotal/Ittotal) *SQ(Isp) *SQ(sigma): 


Qp ~ 2.0*Ecphype([j) *Istotal*(1.0-(Istotal/Ittotal)) 
-hsq* (Istotal/Ilttotal})*(1).0-(lstoetal? it to.11)) 
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+(Istotal/Ittotal) *SQ(Is) *SO(siqm) ; 


arg - ((Ittotal/(Ittotal-Istota))) 
*((Isp*sigma)+ (signg‘sqrl (O))))/Crqut (head). 
if (arg > 0.99999999) arg-1.0: 


etah[j]) © §7.2957795131*acos (arg) ; 
argp = ((Ittotal/ (Ittotal-Istotal)) 
*((-Is*sigma) #(signgp*sqrt (Op) )))/dsqet Cheght 
if (argp > 0.99999999) argp-1.0: 
etahp[j]) = $7.2957795131*acos(argp);: 
Ecexp[)]) - (Ec{lJ)-Ecfinal) *exp(ecfactor*xp(j]) tRefinal: 


zcexpdel [ j) = Ecexp[j)-Ecexp[1]; 


Ecpexp[j]) (Ecp{1)-Ecpfinal) *exp(ecpfactor*»p[j]) tRepfinal: 


Ecpexpdel[j: = Ecpexp[j]-Ecpexp([1)- 


Q = 2.0*Ecexp[j) *Istotal*(1.0-(Istotal/ittotal)) 
-hsq* (Istotal/Ittotal)*(1.0-(Istotal/Ittotal)) 
+(Istotal/Ittotal) *SQ(Isp) *SQ(sigma); 


Qo = 2.0*Ecpexp[j])*Istotal* (1.0-(Istotal/Ittota))) 
-hsg* (Istotal/Ittotal)*(1.0-(Istotal/Ittct4a!)) 
+(Istotal/Ittotal) *SQ(Is) *SQ(sigma) ; 


arg e ((Ittotal/ (Ittotal-Istotal)) 
*((Isp*sigma)+(signg*sqrt (0))))/ (sqrt (haq)): 
if (arg > 0.99999999) arg=1.0; 


eca[j) e 57.2957795131*acos (arg); 


argp e ((Ittotal/(Ittotal-Istotal)) 
*((-Is*sigma)+(signqp*sqrt (Qp))))/ (aqri (ay) ): 
if (argp > 0.99999999) argp=1.0: 


ezap[}j) = 57.2957795131*acos (argp): 


c 


peeerveanaar initialize graphics eeennewaaeea / 


init(); color_scale “cyanblue") ; 

Grey _scale(“gGreyscale-"); window0(); 

bgcol(7); erase(); <zolor(0): 

move (220,-110); prinzi("MINGORI DUAL SPINNER MASS SPRING SYSTEM"): 


window (22,-85,979,872 : bgceol(3): 
erase(); scale(0,1072°35,10000,0); rect (0,0,10000, 10000): 
e1ize(1,1):; 

~ove (20, 110+200); date (); 

~sve(-160,110+70); printf ("0"); 

ove (10,110-90); printf ("w3 - w3 initial”:: 


vector(0, .5,110,19); -sve(-160,15+70); printf£("0"); 
move (10,194+200); printf ("w2"); 


vector(9,15,110,18); -zve(-160,18+70); printf£("0"); 
~ove (10,18+200); printf (“wl"); 
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vector(0,17,110,17): move (-160, 17470): 


pringei( 0"): 


move (10,17+200): print£{("“z prime dot"): 


vector (0,16,110,16): move (-160, 16+70): 


print£ ("0"); 


move(10,16+200): printf£("z prime"); 


vector(0,15,110,15): move (-160, 15+70): 


print£("0"): 


move (10,15+200);: print£("z dot"): 


vector(0,14,110,14): move (-160, 14+70): 


Srint sev): 


move (10,14+200); printf ("2"); 


vector (0,13,110,13): move (-160,13+70); print£("0"): 
move (10,13+200); printf£(“ke-ke 1 (r), 


vector(0,12,110,12); move(-160,12+70); printf("0")-: 


move (10,12+200); printf£("Ec-Ec i (r), 


vector(0,11,110,11): move (-160,11+70); printf("0"); 


move (10,11+200):; printf£("Ecp-Ecp i 


vector(0,0,110,0): 


vector (0, 19+750,11, 19+750); 
vector (0,19+500,11, 19+5S00); 
vector (0, 19+250,11, 19+250);: 
vector (0, 19-250,11, 19-250); 
vector (0,19-500,11, 19-500); 
vector (0,18+250, 11, 18+250);: 
vector(0,18-250,11, 18-250); 
vector (0,18-S00,11,18-500);: 


xscale =~ 10000.0/xp[kount]); 


olor (1)%: 
for (jel: j<kount; j++) 
( 


vector ( (int) (xscale*xp(j)), (int) (13+(0.01l*etotaldel[j))). 


/**e*e* blue - total energy, 


move (-160,0+70); printf£("0"): 
move (10,0+200); printf(“theta, eta (b), etap (r), 


vector (19, 19+750, 110, 19+750): 
vector (19,19+500,110, 19+500); 


vector (19, 19+#250,110,19+250): 


vector (19,19-250,110,19-250): 
vector (19, 19-500, 110,19-500): 
vector (19, 18+250, 110, 18+250): 
vector (19,18-250,110, 18-250): 
vector (19, 18-500, 110, 18-500): 


Echype-Echype i 


(rc), Ecphype-Ecphype i 


rotor hype *****/ 


I aaa Ts 


Ch): 


phen ds 


h- Wh oinit iol”); 


(int) (xscale*xp[jt+lj), (int) (13+(0.0l*etotaldel[j+1)))): 


vector ( (int) (xscale*xp(j]), (int) (13+(0.1*etotaldel[j}))., 


(int) (xscale*xp([j+1J), (int) (13+(0.1l1*etotaldel[(j+1)))): 


vector ((int) (xscale*xp([j)), (int) (12+(50.0*Echypedel[j)})). 


(int) (xscale*xp[j+1]), (int) (12+ (5S0.0*Echypedel {| j+1)))): 


vector ( (int) (xscale*xp[j}), (int) (0+(100.0*%etah[j})), 


(int) (xscale*xp([j+1)), (int) (0+(100.0*%etah[j+1)}))): 
vector ((int) (xscale*xp[j)), (int) (11+(1000.0*etah([j))), 


(int) (xscale*xp([j+1]}), (int) (11+(1000.0*%etah[j+1)}))): 


eEqloz (5); 
for (j=l: j<kount; j++) 
{ 


veccor( (int) (xscale*xp[j}), (int) (11+(50.0*Ecpexpdel[j])), 


/***** purple - platform exponential*****/ 


(int) (xscale*xp[jt+1]), (int) (11+ (50.0*Ecpexpdel {| j#1j))): 


vecior( (int) (xscale*xp[j)), (int) (0+(100.0*etap[j))), 


(int) (xscale*xp[j+1}), (int) (0+(100.0*%etap[jil]}))): 
vector ((int) (xscale*xp[j)}), (int) (11+(1000.0*etap[j})), 


(int) (xscale*xp[j+1]), (int) (11+(1000.0*etap[jtl}))): 
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) 


color (6); jose vellow = rotor exponential *****/ 
for (jel: j<kount: je-! 
( 


vector ( (int! ‘xscale*xp(jJ), (int) (124 (50.0*Ecexpdel[j})). 
(int} (xscale*xp(jt+l)), (int) (12+(50.0*%Ecexpdel (pil yily: 


vector( (int) ‘xscale*xp[j]), (int) (0+(100.0*etal[j))). 
. (int) (xscale*xp({j+1)), (int) (0+(100.0%etal j+1}))): 
vector ( (int! ixscale*xp(j)), (int) (11+(1000.0*etal[j))). 
(inz' (xscale*xp(j+1)), (int) (114 (1000. 0*%etal[jil)))): 


] 


color (0);- /***** black - platform hyperbolic *****/ 
for (jel; j<kount; j+e; 
{ 


vector ((int) (xscale*xp[j)), (int) (11+(50.0*Ecphypedel[(j})). 
(int! (xscale*xp({j+l)), (int) (11+ (50.0*%Ecphypedel(jtlj))): 


vector ((int) (xscale*xp[j]), (int) (0+(100.0*etahp([j))). 

(int: (xscale*xp[j+1lJ]), (int) (0+(100.0%etahp(j+])))): 
vector ((int) (xscale*xp(j]), (int) (11+(1000.0%etahp[j))), 

(inti (xscale*xp{jt+l)), (int) (11+ (1000.0*etahp[j'))))): 


) 
color (4): /***** red = rotor ****e/ 
for (j=-1; }<kount; j++} 

( 


vector ((int) (xscaie*xp[j)), (int) (110+ (1432.39448783* (yp(3) (j)-yrf3y0b)))). 
(int) (xscale*xp($+1)), (int) (110+ (1432.39448783* (yp(3]) [j+l)-ypl3!il))))): 


vector ((int) (xscale*xp{j)), (int) (194+ (1432.39448783*yp(2)(j))). 
(int) (xscale*xp(j+l)), (int) (19+ (1432. 39448783*ypl[7J(jtdi)))- 
vector ((int) (xscale*xp{j)), (int) (18+(1432.39448783*ypl(1)(j))), 
(int) (xscale*xp([j+1)), (int) (184+ (1432.39448783*ypl))titlbii): 


vector ((int) (xscale*xp[j)), (int) (17+(1000*yp[{7)[j))). 
(ant) (xscale*xp[j+1)), (int) (17+ (1000*yp(7)(j41)))): 


vector ( (int) ixscale*xp(j)),. (int) (16+(1000*yp(6){j))) 
(int) ‘xscale*xp(j+1l])), (int) (16+ (1000*yp(6)[j+1}))): 


ve tx ( (int) ixscale*xp(j)), (int) (15+ (1000*yp(5)[(3))). 
(int) ‘xscale*xp[{j+1]), (int) (15+ (1000*yp(5)[(j+#1)))): 


vector ((int) ixscale*xp(jJ)), (int) (14+(1000*yp(4) (jJ))., 
(ant; ‘xscale*xp[j+1l)), (int) (14+ (1000*yp(4)[j'1)))): 


vector ( (int) (xscale*xp[j)), (int) (13+(0.01*kedel[j])), 

(int) .xscale*xp[(jt+l)), (int) (13+ (0.01*kedel[jt+1)))); 
vector ((int) (xscale*xp{j]), (int) (13+(0.1*kedel{j))), 

(int) ‘xscale*xp(j+l)), (int) (13+ (0.1*kedel[j+1]}))): 


vector ( (int) {xscale*xp(j)), (int) (12+(50.0*Ecdel[jJ)). 
(int) .xscale*xp[j+1])), (int) (12+ (50.0*Ecdel[j'l]))): 


vector ((int) (xscale*xp[{j]), (int) (11+(50.0*Ecpdel[j))), 
(int) :xscale*xp[j+1J), (int) (114+(50.0*Ecpdel[j'l))}): 


vecioz( (int) (xscale*xp[{j)), (int) (0+(1000.0*hdel[j))). 
(int) ‘zscale*xp[j+l]), (int) (0+(1000 .0*hde] (j'1]))): 
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vector((int! (xscale*xp{|j)), 
(int) (xscale*xp[j+1 
vector((int) (xscale*xp[j)), 


(int) (0¢(100.0*thetafj))). 
J), Cint) (04 (100.0*thetal jelyyio: 
(int) (114(1000 .0*theta[j))). 


(int) (xscale*xp[j+1]), (int) (114 (1000. 0*thetaf jrldiid: 


/* vector( (int) (xscale‘*xp[j)), 
(int) (xscale*xp(j+1 

o/ 

/* vector ((int) (xscale*xp[j)), 
(ini) (xscale*xp[ j+1 

a] 

/* vector ((int) (xscale*xp[j)), 
(int) (xscale*xp[j+1 

i 


e ) 
window) (): 
color(Q); 


frre 
move (380,875); prints (" 
/* 


move (22,910): 


move (500,910); printf£(“eps = te 


move (22,940): 


print initial, 


INITIAL CONDITIONS 


(int) (04(0.0001*ke[j})). 
1), (int) (0#(0.0001*kef[ jt1)))): 


(int) (0+(0.0001l*etotal[j})). 
J), (int) (0+(0.0001*ctotal{ jrlj))): 


(int) (0+(0.1*%h[j])), 
J). (1NO) (0400s Tshirt): 


anal 


final data 


kount -— Yd", hovweat J: 


printf("time =-%3.0f£ seconds", xp[kount)): 


td/td",eps,nbad, nok): 


print£ ("Wl <04.2£ W2 =84.2£ W3 <04.2£", yp(l] (1), ypl2) (1). yrtIltl)): 


move ($00,940): princft("w3p= %t9.7f 
move (22,980); 
move (22,1010): 
move (22,1040): 
move (22,1070); 


princf£("I1l =84.0f 
printf("M = %4.0f 
printf£("a =84.1f 
printf("z =-84.2£ 


move (500,980); printf("Ilp=84.0f 
move (500,1010); print£("Mp =-44.0f 
move (500,1040); printf ("ap=84.1f 
move (500,1070); printf(“zp =%4.2f 


move (22,1110); 
MOVE(22-21059)> 
move (22,1160): 


printf("h i#t9.4f 
princf(“h f£-t9.4¢£ 
prints (“hdelper=t5. 


L = ¥3.1£",w3p, Lb): 


I2 —84.0£ I3 =84.0£",11,32,1): 
m =-&7.5£",M,m); 

k =~84.1f C ete lia, k. Cc): 
zdot\=84.2f",yp(4)(1l)].yp(5) (1)): 


I2p-t4 .0f Tapet4 Of" Tipp.) ph) 
mp=%7.5f£",Mp,mp) : 
kp=84.1f cp=t4.1f",ap, kp, cp) : 


zdotp = 2f", yp(6)(1),ypl7)11)): 


ke i=87.1£",h[1),ke[1]): 
ke f£-87.1f£",h[kount],ke[kount)): 
2£",100.0* (h[kount)J<-h[1))/hi1)): 


move (250,1160); printf (“kedelper=85.2£",100.0% (ke[kount]-ke[1l]))/hefl)l: 


move (500, 1160): 


printf(“e celper=85.2£",100.0*(etotal[kount])-etotal[1l}) /ototal{))): 
move (730,1160); prinzf£("Ecp delper=%5.2£",100.0* (Ecp{kount)-Ecp{!))/Ecptil)): 


move (§00,1110): 
move ($00, 1135): 


Prinzf("Ec isWS sae 
printi("Ec f£=t9.4f 


theta i=37 .S£",Ec(1),thetal[l]): 
theta f=37.5f£",Ec[kount),thetaftbount ]); 


=/ 
move (22,900); Dine cit - platform roto? +). 
move (22,920): prints ("lip "= 26.267 1p: 


move (240,920); printf£("Il - %6 
move (22,940); printt(“I3o5 = %6. 
move (240,940); printf("I3 = 6. 
move (22,960); prin. (Mp = %6. 
move (240,960): printst("M = %6. 
move (22,980); Princ. (mp = &6. 
move (240,980); printf ("m = %6. 
ove (22,1000)> princs(7ap = %6. 
move (240,1000); printf("a - %6. 
move (22,1020)3" (prine<( ko = %6. 
move (240,1020);: printf ("k = %6. 


LZ Eas) 


Zee LOD) * 
2f°213)2 
2£oMp): 
ZEN 
2£° MD) .2 
Ze am yas 
Zt" apc 
Ztail: 
2£""Ko) 
7.1) Gaara ea 
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move (22,1040): printf£("“cp = eh. 2E eps 
move (240,1040); printf("“c = &6.20°", ¢)": 
move (22,1060); printf("L 6.40",1L); 
move (240,1060); printf("Is/It = $6.3f£",Istotal/Ittotal): 


move (22,1100): 
move (240,1100); 
move (22,1120); 
move (240,1120): 
move (22,1140) ;- 


printf ("zp 
printf ("z 
printf (“dzp 
printf (“dz 
printf (“w3p 


@beeetenadt 
rd 
fon) 


move (240,1140); printf ("“w3 &6 
move(22,1160); printf ("wi %6 
move (240,1160); printf (“w2 %6 


move (480, 900): 


ver, ¥PLoO) (1}) 
-2f",yp[4)(1)) 

met yp tT il) h: 

eet YN LOL) 

.2f", yp(3)[1)+sigma): 
Bet, Yous) (1) }: 

eet eeypli) (1i)): 

mee yp(2)[1))-: 


printf("t3.0f seconds", xp[kount]); 


move(730,900); printf(“h i - 89.4f",h[1)): 
move (480,920): printf("td / td / %td",nbad, nok, kount): 
move (730,920): printf("“h f = 89.4f",h[kount]);: 


move (480,940): printf ("eps 


move (730, 980); 

move (480, 1000): 

move (730, 1000): 
“printf ("edelp 


printf("“e f 
printf (“kedelp 


move (730,940); printf("“hdelp - $5 
move (480,960); printf("ke 4 - 69 
move (730,960): printf("e i - 4&9 
move (480,980): printf("“ke f = 4&9 


move (480,1030); printf ("“Ecpi = $9. 
move (730,1030); printf£("Ec i = $9. 
move (480,1050);: printf ("“Ecpf = 69. 
move (730,1050); print£("Ec f = $9. 
move (480,1070): printf("“lambdap = 8&5. 
move(730,1070);: printf("“lambda = %5. 


move (480, 1090): 
move (730,1090):; 
move (480,1110); 
move (730,1110)-;> 


printf("signgp «= 
printf (“signg 


printf (“dEc/lambda 


move (480,1130): 
move (730,1130); 


printf (“Ecfp 
printf ("Ecf 


move (480,1155): 
move (620,1155S); 
move (760,1155); 
move (480,1175); 
move (620,1175): 


move (760,1175); printf("eta 


{eee 
move (600,-33):;: 
move (600,-13): 
move (600, 27); 
move (600, 47); 


move (980, -33): 
move (980,-13): 
move (980, 7): 
move (980,27); 
move (980, 47): 
move (980,67): 


printf (“td", nbad): 
prints ("td", nok): 
Printf("%d",kount); 
printf (“nbad"): 
printf (“nok"); 
printf ("kount"); 


move ($80,S$7): printf("%6.4£", 


=—to.35f", 100. 


tS. 
No: 
printf ("dEcp/lambdap- NG.3f", 
= $6.3£", (Ec[kount)-Ec[1])/lamrted-): 


= $8. 
= %8. 


upper right hand corner 
printf(“Ec initial/final 
printi(“Ecp initial/final = %g / tg",Ecp[1l),Ecp[kous:j); 
Pprintf£(“Ec factor / Ecf 
printf(“Ecp factor / Ecfp = %tg / %&g",ecpfactor,FEclfp): 


(100.90* 


VW7e",eps); 

Pat" 100.0* (ht kounViehilijattit tl). 
.4f£",ke[1]): 

.4f3% 


etotal[i]}): 
4f",ke[kount])):; 
4f"“,etotal[kount]}): 


3£", 100.0* (ke[kount] -ke[)])/Fel)i): 


O* (etotal[kount)-etotal{]]})) fer et algij). 
Sf, CPL) |); 

Sf EG) i 

3f£",Ecp[kount]): 

3f£",Ec(kount)): 

3£", lambdap) : 

3f£", lambda) ; 

Ze, Signgp) : 

2£", signq): 

(Ecp[kount )-Ecp[1)}) /lamsdap) : 


6£", Ecfp)- 
6f£",Ecf): 


printf(“theta i =%6.2f£",theta[1l}): 
printf("etap i ©%6.2f",etap[1))-: 
printf(“eta i =%86.2f",eta[1l]}): 

printf ("theta f =86.2£",theta[kount])): 
printf("etap £ -%86.2f£",etap[kount));: 

f -“%6.2f",eta[kount]}): 


een / 


=tg / Ag".Ec(1),Ec(kount)): 


e %g / tg",ecfactor,Ecf}: 


(ke[kount)-ke[1)))/he[1l)): 
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move ($80,117); printf ("¥6.4£", (100.0* (etotal [kount J-etotal (I) ))/-rerartly: 
move (980,137): print£("%86.4£", (100.0° (h{kount p-hfl))) ZbGl pd: 


move (980,157): printf ("kedel"): 
move (980,177); printf("“edelp"); 
move (980,197): printf(“hdelp"): 


move (980, 227); print£("tg",theta[1l}); 
move (980,247); printf ("%g",thetal(kount}): 


move (980,267);-printf("th 1"); 
move (980,287); princf("th £"); 


move (980,317): printf (“8g",etafll))-: 
move (980,337): printf ("%tg",etalkount]): 


move (980, 357): prinzfi ("eta 1"); 
move (980,377); princf("eta £"); 


move (980,407); printfi("tg",etap[1)): 
move (980,427): printf (“%tg",etap[kount}): 


move (980,447); printf ("etapl"): 
move (980,467); printf("etapf"); 


move (980,497); print£("%tg", Ec(kount]-Ec[1))-: 
move (980,517); printf ("%tg", lambda); 
move (980,537); printf£(“%g", signgq); 


move (980,557); printf ("Ecdel"); 
move (980,577); printf ("lambda"); 
move (980,597); printf ("signg"): 


move (980,627); printf ("%tg", Ecp(kount)]-Ecp[1l]}): 
move (980,647); printf ("%tg", lambdap) ; 
move (980,667); prints ("%g", signgp);: 


move (980,687); prints ("Ecpdel") ; 
move (980,707): prints ("lambdap") ; 
move (980,727): print£("signgp"); 


move (980,757): printl("&g", Istotal/Ittotal): 
move (980,777); prints ("Is/It"); 


zee_matzrix (yp, 1,h, 2,MAXARRAY) ; 
free_vector (*p, 1,MAXARRAY) ; 
free vector (ystart,=.N); 
free_vector (ke, 1,MAXARRAY); 
free_vectoz (kedel,1,MAXARRAY) ; 
free_vector (etotal,1,MAXARRAY) ; 


free_vector (etotalde.,1,MAXARRAY) ; 


free_vector (h,1,MAX#RRAY); 

free vector (theta,1,¥AXARRAY) ; 
free_vector (hdel,1,™“AXARRAY) ; 
free_vectoz (Ecexp, 1,“AXARRAY) :; 
free_vector (Ecpexp,1,MAXARRAY) ; 
free_vector (etah, 1,M=XARRAY); 

} 


free_vector (Ec, 1, NAXMAFPAT); 
free_vector (Ecdel, 1, MAXNAPHA7) ; 
free_vector (Ecp, 1,MAAARPAT) ; 
free_vector(Ecpdel, 1, M@AXAFRAY); 
free_vector (Echype, 1, MAXAPRAY) : 
free_vector (Echypedel, 1,A“APREAY) ; 
free_vector (Ecphype, 1,NAXARPAY) : 
free_vector (Ecphypedel, 1, iin“APEAY) : 
free_vector (eta, 1,MAXAPRAY) : 
free_vector (etap, 1,MAXARRAY) ; 
free_vector (Ecexpde)], 1, NAZAPPAY): 
free_vector (Ecpexpdel, i, MAZAKRRAY) ; 
free_vector (ctahp, 1, MAAARAY) : 


i 
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